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Chapter  1 

INTRODUCTION 


The  electromagnetic  scattering  from  a  subclass  of  superquadric  surfaces, 
specifically  two-dimensional,  perfectly  conducting  superelliptic  cylinders  is 
treated.  Scattering  may  be  defined  as  the  modification  of  the  electromag¬ 
netic  radiation  fields  due  to  the  presence  of  complex  geometries. 

Computer  modeling  is  an  active  area  of  high-frequency  electromagnetic 
research.  It  concerns  itself  with  the  construction  of  efficient  computer  pro¬ 
grams  to  calculate  values  for  the  various  antenna  and  radar  cross-section 
parameters  of  complex  antennas  and  scatterers.  Some  of  these  computer 
codes  [14], [15]  are  capable  of  modeling  quite  complicated  antennas  and  scat¬ 
terers  such  as  reflectors,  ships,  aircraft,  spacecraft  and  many  other  actual 
structures. 

The  purpose  of  this  report  is  to  extend  the  scope  of  computer  modeling 
codes  to  include  super  quadric  shapes  in  order  to  represent  complex  geome¬ 
tries  by  a  new  class  of  analytic  functions.  Superquadrics  show  great  promise 
of  providing  researchers  and  engineers  with  a  powerful  family  of  parametric 
shapes  for  geometrical  modeling.  One  of  the  main  problem  areas  today  in 
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modeling  is  the  lack  of  a  unified  mathematical  formalism  and  specification 
language  for  geometrical  objects.  By  providing  modelers  with  the  fprility  to 
generate  a  wide  variety  of  shapes  from  a  small  number  of  intuitive  parame¬ 
ters,  superquadrics  may  prove  to  be  a  step  in  the  right  direction  toward  the 
needed  mathematical  basis.  Since  superquadrics  allow  complex  surfaces  to 
be  generated  and  modified  easily  and  interactively,  there  is  also  hope  that 
they  would  integrate  naturally  with  an  evolving  specification  language  for 
geometrical  objects. 

This  report,  then,  represents  a  small  step  in  examining  the  potential  of 
these  surfaces  for  this  purpose. 
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Chapter  2 

THEORETICAL 

BACKGROUND 


2.1  Introduction 

The  solutions  of  electromagnetic  problems  consist  of  solutions  to  Maxwell’s 
equations  and  the  equation  of  continuity,  together  with  appropriate  bound¬ 
ary  conditions.  There  are  three  high-frequency  techniques  of  particular 
interest  here,  Physical  Optics  (PO),  Geometrical  Optics  (GO)  and  the  Uni¬ 
form  (Geometrical)  Theory  of  Diffraction  (UTD). 

2.2  Geometrical  Optics  (GO) 

Geometrical  Optics  is  an  approximate  technique  tha'  can  be  used  to  rep¬ 
resent  radiated,  reflected,  and  refracted  fields.  GO  can  be  derived  via  an 
asymptotic  series  (Luneburg-Kiine)  solution  of  the  Maxwell’s  equations;  the 
leading  term  of  the  series  is  the  GO  field. 

The  GO  field  is  discontinuous  across  a  shadow  boundary.  Its  amplitude 
is  governed  by  the  conservation  of  energy  in  a  ray  tube  as  it  travels  along 
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the  ray  path.  The  phase  is  proportional  to  the  length  of  the  ray  path,  and 
ray  tubes  are  defined  by  surfaces  normal  to  the  ray  path  through  which  the 
flow  of  power  is  a  constant.  GO  fails  when  the  energy  of  a  ray  tube  must 
pass  through  a  point  or  line.  Such  points  and  lines  are  called  caustics,  aud 
they  signal  the  attempt  by  GO  to  represent  the  flow  of  a  finite  amount  of 
power  through  a  vanishing  area.  In  two  dimensions,  the  GO  reflected  field 
is 


El(rl  =  ~E\{Qt)  •  ^ 

/ .  PT  c-ik.' 

(2.1) 

l(PT  +  sr) 

Kif)  =  +Hi(Qr)  ■  ( 

1  pr  c~jh‘r 

(2.21 

\l(pr  +  3')C 

where 


Qr  =  reflection  point 
E\(Qr)  —  incident  Et  field  at  Qr 
H\(Qr)  =  incident  Ht  field  at  Qr 

s'  —  distance  from  Qr  to  source 

ar  —  distance  from  Qr  to  receiver 

pT  =  caustic  distance  for  the  reflected  ray 

pT  p'  +  Rc(Qr)  cos  0* 

—r  =  caustic  distance  for  the  incident  ray 
P‘ 

8 ‘  =  angle  of  incidence  =  cos_1(— J*  •  n) 

=  principal  radius  of  surface  curvature  at  Qr. 
h  =  normal  to  the  surface  at  Qr 
t  =  tangent  to  the  surface  at  Qr 
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i’  =  direction  of  the  incident  ray 
aT  —  direction  of  the  reflected  ray 
r  —  direction  of  the  source 
r  =  direction  of  the  receiver 

as  can  be  seen  in  Figure  2.1.  The  direction  of  the  reflected  ray  is  defined 


Figure  2.1:  Geometry  for  the  GO  Reflected  Field 
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by  the  law  of  reflection  and  is 


n-iT  =  -*■*<.  (2.3) 

The  point  of  reflection,  Qr,  is  a  point  on  the  surface  such  that  the  law  of 
reflection  is  satisfied.  For  far  field  scattering  from  Qr,  h  ■  f  —  h  ■  r' . 


2.3  The  Physical  Optics  Procedure  (PO) 


The  currents  induced  on  a  general  scatterer  are  unknown.  If  the  true  cur¬ 
rents  were  known,  then  exact  field  values  could  be  calculated  using  the 
radiation  integrals  [9j.  Physical  Optics  is  a  procedure  where  unknown  cur¬ 
rents  are  approximated  by  equivalent  currents  based  on  the  incident  G.O. 
fields.  The  equivalence  theorem  allows  replacement  of  the  original  scatter¬ 
ing  geometry  and  the  actual  surface  currents  by  equivalent  surface  currents 
flowing  in  free  space.  These  approximate  G.O.  currents  are  then  used  to 
calculate  the  scattered  fields.  The  currents  induced  on  the  surface  of  a 
perfect  electric  conductor  arc  assumed  to  be 


f  f  2n  x  U\ 

'  "  t  0,  i 


in  the  Lit  region 
in  the  Shadow  region 


(2.4) 


where  n  is  the  unit  normal  to  the  surface. 

Physical  Optics  is  useful  because  the  form  of  the  assumed  currents  is 
simple  and  the  resulting  integrals  lend  themselves  to  either  numerical  or 
asymptotic  (high-frequency)  analysis.  It  is  not  possible  to  rigorously  state 
the  conditions  of  validity  of  PO,  but  several  guidelines  may  be  employed 
when  determining  its  validity  in  a  given  problem.  In  general,  whenever  the 
assumed  currents  are  not  a  good  approximation  to  the  actual  currents  then 
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PO  will  be  invalid  somewhere.  For  example,  in  a  bistatic  scattering  con¬ 
figuration  where  currents  in  a  shadow  region  contribute  significantly  to  the 
total  field,  the  physical  optics  approximation  is  no  longer  accurate.  Another 
source  of  error  in  the  PO  formulation  is  the  assumed  termination  of  G.O. 
currents  on  the  scatterer  surface.  Since  the  actual  equivalent  currents  do 
not  end  abruptly,  evaluation  of  the  PO  integral  yields  contributions  from 
the  endpoints  of  integration  which  are  nonphysical.  When  these  spurious 
current  termination  contributions  can  be  identified  and  removed,  the  PO 
result  becomes  more  accurate.  Care  must  be  exercised  however  when  de¬ 
termining  whether  a  term  is  a  false  current  termination  or  the  result  of  a 
bona-fide  discontinuity  with  physical  causes. 

When  PO  is  a  good  approximation,  and  if  the  stationary  phase  condition 
is  applicable,  then  a  recovery  of  the  GO  result  is  possible  from  PO.  Because 
PO  is  a  spatial  integration  of  surface  fields,  it  produces  bounded  results  in 
situations  where  the  conditions  required  for  a  valid  GO  result  do  not  hold. 
This  ability  of  PO  to  treat  scattering  not  possible  with  GO  suggests  that  it 
represents  a  viable  and  useful  format  for  many  problems.  This  is  the  main 
reason  that  PO  was  chosen  to  characterize  the  scattering  from  superquadric 
surfaces. 


2.4  Uniform  Theory  of  Diffraction  (UTD) 

The  Uniform  Theory  of  Diffraction  [10]  is  a  uniform  version  of  the  Geomet¬ 
rical  Theory  of  Diffraction.  The  GTD  is  an  extension  of  Geometrical  Optics 
that  postulates  the  existence  of  “diffracted  rays”.  Recall  that  the  GO  field  is 
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discontinuous.  Since  actual  fields  must  be  continuous,  then  diffracted  fields 
are  generated  that  eliminate  this  discontinuity.  These  diffracted  fields  are 
added  to  the  GO  fields,  i.e.,  Uft4  =  U^,  4-  t/aif. -  The  postulates  of  GTD  are 
very  similar  to  those  of  GO. 

1)  The  ray  paths  may  be  found  as  a  generalization  of  Fermat’s  principle- 
diffraction  points  occur  at  places  such  that  the  total  ray  path  is  an 
extremum. 

2)  Diffraction  like  reflection  and  transmission  is  a  local  phenomenon  at 
high  frequencies. 

3)  For  a  diffracted  ray,  power  is  conserved  in  a  tube  of  rays  and  the  phase 
of  the  diffracted  field  is  proportional  to  the  length  of  the  traversed 
ray. 

Generally,  a  UTD  solution  is: 

•  accurate 

•  valid  at  reflection  and  shadow  boundaries  and  in  the  shadow  region 

•  valid  for  an  incident  ray  optical  field  with  arbitrary  wavefront  curva¬ 
ture 

•  computationally  efficient. 
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Chapter  3 

SUPERQUADRIC 
SURFACES  AND  THEIR 
PROPERTIES 


Superquadric  surfaces  are  generalizations  of  the  quadric  surfaces.  The  def¬ 
inition  of  a  quadric  surface  is  the  locus  of  all  points  (x,y,  z)  that  satisfy  the 
equation 

Ax3  4  By 3  4  Cz3  4  Dxy  4  Exz  4-  Fyz  4  Gx  4  Hy  4  Jz  4  K  =  0  (3.1) 

for  arbitrary  constants  A ,  B,  C,  D,  E,  F,  G,  H,  J  and  K.  In  two  dimensions, 
Equation  3.1  is 

Ax3  -f-  By 3  4  Dxy  4  Gx  4  Hy  4  K  =  0.  (3.2) 

Examples  of  quadric  surfaces  include  the  ellipsoid,  the  hyperboloid  of  one 
sheet,  the  hyperboloid  of  two  sheets,  the  elliptic  paraboloid,  the  hyperbolic 
paraboloid,  the  elliptic  cone,  and  the  quadric  cylinder.  Superquadrics  do 
not  satisfy  Equations  3.1  and  3.2  except  for  certain  special  cases.  For  our 
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purpose,  we  restrict  our  attention  to  closed  supereiliptic  and  superellip- 
soidal  surfaces. 

The  transition  from  quadrics  to  superquadrics  is  accomplished  by  al¬ 
lowing  the  exponential  powers  of  x,  y  and  z  in  Equation  3.1  to  take  on 
arbitrary  values.  The  difficulty  of  raising  negative  numbers  to  fractional 
powers  is  resolved  through  the  careful  application  of  absolute  value  signs 
and  sign  functions.  The  proper  combination  of  absolute  value  signs  and 
sign  functions  ensures  that  a  closed  quadric  surface  generalizes  to  a  closed 
superquadric. 

To  illustrate,  consider  the  equation  of  a  quadric  ellipsoid  in  rectangular 
coordinates  that  meets  the  x,  y  and  z  axes  at  ±o,  ±6  and  ±c. 

(!)'+(?)'+(!)’=«• 

It  is  useful  to  write  the  surface  equations  in  parametric  form  so  that  u  and 
v  are  the  principal  directions  of  the  surface.  The  position  vector  is  then 

f  =  xa  ■  sinu  cosv  +  yb  •  sinu  sin  v  +  zc  •  costt  (3-4) 

where 

0  <  u  <  7r  and  0  <  v  <  2 ff.  (3.5) 

A  normal  to  the  surface  is  given  by 

_  *  y  z 

n  =  —  •  sin u  cos  v  +  —  •  sin  u  sin  v  +  -  •  cos  u.  (3.6) 

a  b  c 

Generalizing  the  quadratic  exponent  in  Equation  3.3  and  taking  note  of 
the  proper  absolute  values,  we  obtain  the  equation  for  a  superellipsoid  in 
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rectangular  coordinates, 


'  X 

y 

I'll  Vltv\ 

Z 

+ 

+ 

— 

.  a 

b 

c 

for  which  the  components  in  parametric  form  are 

x  =  a  •  |sinup^  •  |cos  *  sign  [sinti  cost;] 

y  =  b- |sinuj*2/^  •  |sint;|^,/1^  •  sign  [sin u  sinw] 

z  =  c  •  jcosup^  •  sign  [costt] .  (3.8) 


Note  that  in  three  dimensions  the  superellipsoid  can  have  two  different 
“squareness”  parameters  vx  and  t/j  applied  to  the  principal  directions  u 
and  v.  Note  also  that  the  sign]]  function  restores  the  proper  sign  removed 
by  the  absolute  value  operators.  For  brevity,  the  use  of  absolute  values  is 
dropped  in  presenting  the  normal  to  the  surface  of  the  superellipsoid, 

n  =  — (8inu)2~2/,l^(co8t;)a_,/,l,,  +  ~(sinu)2 "2^(sm»)2 ~2^  +  -(cosu)2-2^. 
a  b  c 

(3-9) 

Equations  3.8  and  3.9  along  with  similar  expressions  for  superhyperboloids 
are  presented  in  [2]. 

A  second  parametric  form  for  the  surface  vector  f(u,u)  which  satisfies 
a  quasi-superelliptic  equation 


X 

Vl 

y 

Z 

a 

+ 

b 

+ 

C 

=  1 


(3.10) 


similar  to  Equation  3.7  is  given  by 

.  o  •  cos  ^  (sin 


[((cos^f  +  (sinV»r)((cos  ip)”  +(sinv?)KI)] 


vaviiM 


+ 
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+  y 
+  i 


b  •  sinijj  (sin^)^1^ 


[((cos^)*"1  +  (sin  ((cosy?)1*  +  (sinv?)*^)]1^1 
c  •  cos  if 


+ 


(3.11) 


[(cosy>)^  +  (sinvj)^]1^ 
where  <p  and  ij)  are  used  in  place  of  u  and  v  to  indicate  the  alternate  pa¬ 
rameterization.  This  equation  does  not  generate  the  same  class  of  surfaces 
as  Equation  3.7,  but  it  is  presented  here  as  yet  another  generalization  of 
Equation  3.1  with  possible  value  in  surface  modelling.  For  two-dimensional 
curves  generated  in  the  x-y  plane,  they  are  equivalent. 

3.1  The  2-D  Superquadric  Cylinder 

In  this  section  the  equations  related  to  the  superquadric  cylinder  are  given. 
The  two-dimensional  equation  for  the  surface  of  the  superquadric  cylinder 
shown  in  Figure  3.1  is 

irr 


/(*,!')  = 


+ 


=  1. 


(3.12) 


As  a  parametric  function  of  t  where  —1  <  t  <  +1,  the  surface  may  be 
written  as 

F(<)  =  X(t)x  -|-  Y(t)y 

where 

X(t)  =  ±«(1  -  |<n1/*',  Y{t)  =  bi. 

The  normal  and  tangent  vectors  to  such  a  surface  are 


(3.13) 


(3.14) 


V/  iV  ffl1-11  ■  sign  (3f)  4-  \Y  I*-11  •  .ign  (Y) 

n{x'Y)~wrr  () 
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i(X,Y)  = 
so  that 


and 


Figure  3.1:  Various  superellipses  of  v  —  2,3,10. 

*  -  ~iaV  iri°'~1)  •  si8n  (y) + w  •  si«n  (*) 

yjb”'  |A'|a(‘,-1)  +  oJ*'|r|J(,/"1) 

f(<)  =  ±*«(1  -  IC)1^  +  fa 

n(t)  =  +  (j) 

^{1  -  |<n1(1"1/,')  +  a* 


(3.16) 

(3.17) 

(3.18) 
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Using  the  first  angular  parameterization, 

r  =  xa  |cosi>|^2/^  •  sign  (cos  v)  +  yb  |sint;|^2^  •  sign  (sin  u)  (3.19) 
,  6*  |coswj(2_2/l/)  •  sign  (cost?)  4-  ay  [sin  v|(2-2/,/)  ■  sign  (sin  v)  ,0 

«(V)  — — - »mirais,r;,aa:.-.m 1 .1  -v:;1..,:;1  .  ss - .  10.2U) 

Note  that  the  parameter  v  is  not  the  same  as  the  angle  <f>  seen  in  Figure 
3.1.  The  connection  between  the  angle  <f>  and  parameter  v  is 


v  =  arctan 


[G-*r 


<f>  =  arctan  -  (tan  i>)^ 
a 


(3.21) 


(3.22) 


In  the  second  angular  parameterization,  the  relation  between  and  <f> 


where 


r  a 

V’  =  arctan  j^-  tan  <f> 


<t>  =  arctan  -  tan  ^ 
a 


xa  cosip  -f  yb  sin  ip 

(icosv>r  +  i8in^ti,)1<,l/ 


(3.23) 


(3.24) 


(3.25) 


*6|cos^>|1/  1  •  sign  (cos  ^)  +  ya  jsin  ^\u  1  -sign  (sin  V>) 
t/^jcos^i2^-  ^  +  o2  jsin^|2^-1^ 


(3.26) 


For  a  plane  parametric  curve  of  the  form  of  Equation  3.13,  the  radius 
of  curvature  is  given  by 


[(*'(()]’  +  irMlf'’ 

|X'(t)y»(t)-F'(e)X"(t)r 


(3.27) 
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Thus  the  radius  of  curvature  for  the  superquadric  cylinder  is 


( a3v  sin  <j>7v~ 2  +  b3v  cos  ») 


iW)  = 


(v  —  1)  (o6)(*'_1)(co8^sin^)(,'~J)((«sin^)*'  +  (6cos^)l/)(1+1/*') 

(3.28) 


and 


(aa  sin  il)2v~2  +  b 7  cos  rj)2v  a)WJ) 


W)  = 


ab(  v  —  1)  (costf>anif}Yv~7^(smi/)t/  +  cos  ^*')(1+1/*/) 


(3.29) 


Note  that  if  v  equals  anything  other  than  exactly  2,  the  denominator  of 
Equation  3.29  vanishes  at  <j>  =  nir/2 ,  n  =  0,  ±1,±2.  Equation  3.29  is 
plotted  in  Figure  3.2.  This  is  equivalent  to  a  local  zero  in  the  curvature  at 
the  poles  of  the  cylinder.  The  zero  in  curvature  produces  a  “cusp”  behavior 
at  the  poles  seen  in  Figure  3.3.  This  cusp  in  curvature  is  the  source  of  the 
singularity  when  using  GO  to  calculate  the  reflected  fields. 
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Figure  3.2:  Superellipse  radius  of  curvature  Re(il>)  with  a  =  2  and  6  =  1. 
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Figure  3.3:  Superellipse  curvature  cusp  with  for  v  —  2.0, 2.1. 
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Chapter  4 

GO  AND  PO  SCATTERED 
FIELDS 


This  chapter  presents  the  GO  and  PO  formulations  for  the  scattered  fields 
from  a  superelliptic  cylinder.  By  applying  UTD  concepts,  the  total  scatter¬ 
ing  from  the  superellipse  is  seen  to  be  composed  of  several  mechanisms.  As 
a  minimum  these  include  a  reflected  ray  and  a  creeping  wave  around  the 
back.  Since  the  mechanism  of  interest  here  is  the  reflected  ray,  the  other 
possible  mechanisms  are  left  as  suggestions  for  further  investigation. 

4.1  The  GO  Reflected  Field 

For  the  superellipse  geometry  shown  in  Figure  2.1,  the  relevant  parameters 
are 


n  -  ab  If  cos  +  V  8>n 
^  (lasing)*'  -|-  \bcoB<f>\1')1^1' 
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(4.2) 

(4.3) 


Him  = 

= 


pc  = 


n  = 


/5  = 


i’ 


a  = 


H0e~ih,i  (4.4) 

(q&)<1-t'V1'  sin  t2*"2  -f  b2lf  cos  <j>2l'-2)W2) 

(*/  —  1)  (cos  <£sin  <^)(*'-a)((asin  <f>)v  +  (6cos  <£)*')(1+1/*') 


^H(Qr)  COSO1 

xbv  jcos  •  sign  (cos  (f>)  4-  yov  jsin  <j>\' 
\fb2v\coB  +  a2u  |sin  ^>| 


1  •  sign  (sin 
»(*'-!) 


(4.6) 

^4.7) 


ab  cos  9* 


(jasin<^|‘/  +  jfccos  ^|l')1/,,/ 
ab  cos  9X 

(lasin^l"  -f  licos^j1')1^ 
x  cos  9'  +  y  sin  ff 


P~ 


x  cos  9  -f  y  sin  9 

~p' 

P 


(4.8) 

(4.9) 

(4.10) 

(4.11) 

(4.12) 

(4.13) 


The  far-field  GO  solution  for  ff^(r)  is  therefore 

e-jk(r'+r)  jkJcot9i 


Hrz(i 0  =  Hoy/Fc- 


lL 


(i..i,#i«’+i»co. 


(4.14) 


with  <f>,  Rc(4>)  and  9*  given  by  Equations  4.1,  4.5  and  2.3.  As  is  clear  from 
the  form  of  R,.  (Equation  4.5,  Figure  3.2),  the  solution  produces  infinite 
fields  at  the  zero-curvature  points  of  the  superellipse.  This  iB  clearly  a 
failure  of  GO.  The  full  nature  and  reason  for  the  failure  is  not  clear  from 
Equation  4.14,  but  becomes  evident  from  Figure  3.2  when  ^  =  (0,  n7r/2). 
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4.2  The  PO  Scattered  Field 


For  an  arbitrary  2-D  surface,  the  PO  formulation  for  the  TE  scattered  field 
H'  is 

<415) 

where 

p  =  Distance  from  the  origin  to  the  field  point 
r  =  xx  +yy 
a  —  x  cos  9  +  y  sin  9 
J,(r')  =  2n  x  /T(f) 

H'(r)  =  zeih^  =  £eM* '<>•*' +v*™9') 

di  =  Line  integration  dement, 

as  shown  in  Figure  4.1.  Let  the  variable  of  integration  be  the  parameter  t. 
Then 


rV)  = 

xa(  1  -  \t\v)X^v  +  ybi 

(4.16) 

II  il 

s 

o(l  —  \t\v)x/v  cos  9  -f  bt  sin  9 

-xa  sign  (<)  +  yb(  1  -  |^‘/)(1"rl/,/) 

(4.17) 

(4.18) 

vV(  i  -  \tnH^lM + a*  i<i,(p"l) 

J,(r)  x  a  = 

2Hl(r)(h  x  z)  x  a  —  —2 H'x(r)(t  x  i) 

(4.19) 

— 

2ie^“C0*s  +v,inS  ^  x 

(4.20) 

a  sin  9  sign  t  +  bcoa  5(1  — 

(4.21) 

dt  = 

v/j»(i  -  i«n*(i'i/,')  +  o3jti3("_i) 

\ja*  -f  b2  (1  -  |tn2(l_,M  , 

(4.22) 

(i  -  \t\y~^  dL 
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Figure  4.1:  Geometry  for  the  PO  Scattered  Field  from  the  Superellipse 


The  PO  integral  for  the  scattered  H*  field  is  then 

rr.  Fjk  e~>kp  r+l  r  sin0|fj,/~1  signf 
x  ~  V27T  ^  y_i  [a(i-|tn(1"1M 


+  b  cos  6  x 


j^.Co»S')+bt(»in9+»lnS'))  ^ 


(4.23) 


which  may  also  be  written  as 


u-  =  Me—  fl  l 

*  \2ir  y/p  Jo  | 
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(4.24) 
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Chapter  5 

HIGH  FREQUENCY 
ASYMPTOTIC 
EVALUATION  OF  THE  PO 
INTEGRAL 


5.1  Introduction 

The  Method  of  Steepest  Descents  is  a  general  procedure  for  obtaining  ap¬ 
proximations  to  integrals  with  a  large  parameter  k  of  the  form 

I(k)  =  J_f{z)ek«M)dz.  (5.1) 

The  key  to  the  approximation  is  that  significant  contributions  to  the  inte¬ 
gral  will  arise  only  trom  those  parts  of  the  path  P  that  are  local  maxima 
(saddle  points)  of  Re{g(«)}  and  the  endpoints  of  P.  Contributions  from 
the  rest  of  the  path  will  be  exponentially  smaller  and  may  be  neglected. 
The  theory  of  steepest  descent  analysis  is  well  developed  and  fully  treated 
in  [3,6]. 
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A  related  method,  the  Method  of  Stationary  Phase  is  applicable  to 
integrals  of  the  form 

/(*)=  H  f(x)ejhg{a)  dx  (5.2) 

where  again  k  is  the  large  parameter.  The  method  of  stationary  phase  is 
based  on  the  principle  that  rapidly  fluctuating  oscillatory  functions  tend  to 
cancel  under  integration.  Significant  contributions  to  the  integral  will  arise 
where  there  is  a  “stationary”  point  or  a  local  cessation  of  oscillation.  Note 
that  while  Equation  5.2  is  a  special  case  of  Equation  5.1,  the  Method  of 
Stationary  Phase  is  not  a  special  case  of  the  Method  of  Steepest  Descents. 
This  is  true  because  the  contours  of  integration  and  the  analyticity  require¬ 
ments  are  not  the  same  for  the  two  methods.  Both  methods  yield  the  same 
result  for  the  leading  asymptotic  term  if  the  contour  of  one  is  continuously 
deformable  into  the  contour  of  the  other  without  encountering  any  singu¬ 
larities  or  branch  points.  A  good  discussion  of  the  stationary  phase  method 
is  found  in  [3,11]. 

5.2  Topology  of  the  PO  Integral 

Although  Equation  5.3  is  a  stationary  phase  integral,  it  will  be  evaluated 
by  the  method  of  steepest  descents.  The  first  step  of  an  SDP  analysis 
is  to  examine  the  structure  of  the  integrand  in  the  z-plane.  Writing  the 
stationary  phase  integral  4.24  in  the  form  of  Equation  5.1, 

=  <5-3> 
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where 


I(k)  =  h{k)  +  I3(k)  =  f  ft(z)^^dz  +  t  f7(z)ek^dz  (5.4) 

Jo  Jo 

and 

2 

(1  -  zt'f'Z )  +bcosd  (5-5) 

2  V'"0 

+6COS*  (5‘6) 
q,(z)  =  j  (a(l  —  zv)x^v  (cos  0  +  cos 9')  —  bz  (sin#  4-  sin  0')}  (5.7) 

q7(z)  ~  j  (a(l  —  zv)1^(cos9  +  cos0')  4-  5z(sin0  +  sin#*)).  (5.8) 

The  structure  of  the  integrands  in  Equation  5.4  is  shown  in  Figure  5.1.  The 
original  contour  of  integration  is  along  the  positive  real  axis  from  0  to  1. 
The  functions  g,  ,(z)  have  branch  points  at  the  v  roots  of  unity  located  at 

n  =  0,  ±1,  ±2  •  •  •  i/.  (5.9) 

The  associated  branch  cuts  follow  the  contours  of  Re{g(z)}  =  0  so  that 
subsequent  contour  deformations  remain  on  the  same  Riemann  sheet.  In 
order  to  select  the  proper  canonical  integral,  the  behavior  of  q(z)  near  the 
origin  is  of  interest.  The  origin,  in  the  z-plane,  is  where  the  multiple  saddle 
points  are  situated.  In  the  neighborhood  of  the  origin, 

ql3(z)  =  j  ^a(cos0  +  cos#')  ^1 - ^  T  bz  (am  9  +  sin0')^.  (5.10) 

For  both  A  and  /a,  there  will  be  (u  —  1)  saddle  points  (z,u  for  I\  and  z,7i 
for  /2,  j=l,2,3  •••!/  —  1)  that  coalesce  on  the  endpoint  of  integration  at 
the  origin.  This  is  the  primary  contribution  to  the  integral.  The  physical 


fiiz)  —  —  a  sin# 
/j(z)  =  -fusing 


Figure  5.1:  as-plane  structure  of  q  for  v  =  6. 


meaning  of  the  endpoint  of  P  at  z  =  1  is  a  contribution  due  to  the  false 
termination  currents  of  the  PO  approximation,  and  its  effect  on  I(k)  will 
be  ignored  by  deforming  the  original  contour  away  from  this  point  into  the 
valley  region  shown  in  Figure  5.1. 
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5.3  Backseat tered  Field  from  the  Curvature 
Cusp  (the  Pole) 

For  the  special  case  of  backseat  ter  from  the  pole  (0  =  0'  =  0),  Equations 
5.3  through  5.8  simplify  to 

I{k)  =  2J1(ife)  =  2/,(fc)  =  2  ff0{z)ek^dz  (5.11) 

Jo 

f0(z)  =  b  (5.12) 

q„(z)  =  j2a(l  —  zv)l^v.  (5.13) 

It  is  desired  to  approximate  the  structure  of  5.13  near  the  saddle  point 
z,  =  0  with  a  simplified  exponential  structure  which  can  be  integrated  in 
closed  form.  In  the  neighborhood  of  the  origin, 

,,(*)=?  (5.14) 

For  a  transformation  scheme,  we  have 


m 

=  2  I'  f0(z)ekq°{t)dz  = 

Jo 

:  2b  fl  dz 

Jo 

(5.15) 

~  2  r  G0(s)ekT°<l)  ds 

Jo 

=  2 beik,a  ^  e~fc*‘'d$5.  16) 

Jo  ds  .= 0 

where 

=  f<4. 

(5.17) 

=  To(*)  =  = 

:  j2a  -  3V. 

(5.18) 

Equation  5.18  defines  the  transformation  from  the  z-plane  to  the  s-plane. 
The  deformation  of  the  contour  from  the  upper  limit  of  1  in  the  z-plane 
to  oo  in  the  s-plane  is  justified  on  physical  grounds.  It  causes  the  PO 
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termination  currents  to  be  excluded  from  the  approximation.  Expanding 
Equation  5.18  in  a  power  series  about  z,  —  0, 


dz  | 


=  -2ja(v-l)\ 


,'lu 


jins 

=  eJ  * 


A 


..  1 1/*/ 


2  ja 


(5.19) 

(5.20) 

(5.21) 


where  the  choice  of  n  =  0, 1, 2  •  •  •  (1/  —  1)  is  determined  by  the  path  leading 
away  from  zt.  From  Figure  5.1  the  appropriate  choice  for  the  arg  of  the 
steepest  descent  path  is 


[dz 

argU 

hence  n  =  0.  The  canonical  integral  used  in  Equation  5.15  is  given  in  terms 
of  the  gamma  function  and  is 


)  =  -*- 
.J  tv 


(5.22) 


(5.23) 

The  resulting  expression  for  I(k)  is 

™  ~  ^[2;J  /W*«-»r(±)± 

(5.24) 

- 

(5.25) 

The  reflected  field  from  the  pole  of  a  superellipse  (0  =  9‘ 

Equation  5.3  becomes 

=  0)  given  by 

a-  =  s[*2b rfMI  "  P 
*  V  2ir  v  \t/J  I2ka\  yfp 

(5.26) 
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When  v  —  2,  the  supereilipse  becomes  a  regular  ellipse  and  Equation 
5.26  reduces  to 

#;  = 

which  is  the  well-known  result  for  the  backseat tered  field  from  a  2-D  elliptic 
cylinder.  When  v  — ►  oo,  the  supereilipse  becomes  a  box  centered  at  the 
origin  with  width  2b  and  depth  2o.  In  this  case,  Equation  5.26  reduces  to 

H-=2bMe-^e^  (s-28) 

which  agrees  with  the  PO  result  for  the  broadside  backscatter  from  a  strip 
of  width  2b  displaced  along  the  x-axis  by  a.  While  the  analysis  is  carried 
out  for  integer  values  of  i/,  comparison  with  numerical  results  shows  that 
Equation  5.26  remains  valid  for  all  real  values  of  v  between  2  and  oo. 


'b*  e~ikp 
2  a  Jp 


tjk3a 


(5.27) 
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5.4  Reflected  Field  Near  and  Far  from  the 
Pole  for  v  —  3 


The  analysis  for  the  case  when  v  —  3  is  based  on  the  asymptotic  expansion 
found  in  [7]  which  describes  the  arbitrary  configuration  of  two  simple  saddle 
points  situated  near  an  integration  end  point.  To  find  the  saddle  points  we 
set  the  derivative  of  Equation  5.10  equal  to  zero.  Then, 

=  3  ^-a(cos0  +  cos#') 
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- jrrr  T  b  (sin  9  -f  sin  9' 

(1  -  z3)2/3  V 


sin0')^  =  0 

(5.29) 

(5.30) 

(5.31) 


or  for  J1} 

zt  ~|2_  6 (sin#  4-  sin#') 

(1  —  z,3)1^3  a(co80  +  cos0') 

and  for  I2l 

|2  6 (sin 9  4-  sin#') 

(1  _  z,3)1/3  a  (cos  9  4-  cos  6') 

so  that  the  saddle  points  for  It  are  located  at 

r  hn(^)r  r 

~  .  [  IttZ]}*'2  ft,  •  f£±snl3/!  *5'32* 

b[ac°8l  jj]  -liMnlz  j] 

r  hin(^)l3,J  l‘/s 

2”  +  [j  [aco,  (*f  ))3/’  +  [»«n  (“)1S/!J  5'M 

and  those  for  I3  are  located  at 

r  r  / _ _  \  1 3/2  T1/3 


(5.33) 


Itmtqz))*  |‘/3 

r 

[«co«  rfcttii3/a  +  it.in 


(5.34) 
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The  asymmetry  in  the  saddle  point  locations  would  cause  great  complica¬ 
tions  in  the  analysis.  For  relatively  small  values  of  (0  -f  0*)/2,  however,  the 
saddle  points  may  be  approximated  by 


2,  = 


I  3/2 


1/3 


(S±£)]3/>+[krin(S±£)]3'3 

(5.36) 

=  +3  \Z‘\ 

(5.37) 

=  -j  \*A 

(5.38) 

Z2i  =  ~\Z'\ 

(5.39) 

zn  ~  +  \z»  1  • 

(5.40) 

For  larger  0,8',  the  saddle  points  are  widely  separated  and  the  inaccuracy 
of  the  above  assumption  diminishes  exponentially.  Because  of  symmetry  of 
Equations  5.37  through  5.40,  the  following  relations  hold  true. 


fi(zn)  —  /l(zu)  —  /j(2ai)  —  /j(z2j)  (5-41) 

=  b  cos  0  +  a  sin  8  _  ^(i/3))  (542) 

?"(* 21)  =  -9"(*2a)  =  jq"{*n)  - (5.43) 
=  2ja(cos 8  +  cos0')  (^  _  *,)(»/»))  '  (544) 

The  integrals  I\  and  Jj  are  treated  separately,  though  the  strategy  of  anal¬ 
ysis  is  the  same  for  both.  If 


Ir(k)  =  [°°T  fT(z)ehqr^  dz  =  /°°r  Gr{s)ekTr(t)  da  r  =  1,2  (5.45) 

J* r«  J»r* 

where  the  functions  qr(z)  have  two  first-order  saddle  points  at  zri  rJ  then 
the  asymptotic  expansion  of  Jr(fc)  valid  uniformly  as  zri  —*  zri  —*  zra  and 


31 


as  k  — *  oo  is  given  [7]  by 

pfcorO 

/,(*)  ~  ir(/T.(«.i)  +  «,(«.>)]  -j^Gri(>l,k,l\k''3,r.) 

+  *  [#,<»„)  -  ».(»„))  *'/3»™) 


o^Tf(*ra) 


r  =  1,2 


M*\  -  -?«) 

where  the  transformation  from  the  z-plane  to  the  s-plane  is 

a3 

qr(z)  =  rr(s )  -  ar0  +  srjT  — 

qr(zra)  =  Tr(sra) 

«,0  =  J[9r(2l)  +  9r(^2)3 


Vv~r  =  {qr(z l)  -  gr(z2))]  2  Srl  =  ~ar2 


(5.46) 

(5.47) 

(5.48) 

(5.49) 

(5.50) 


H,(>)  = 

«*>z 

dz 

T2srJ 

^  *r\ t*r7 

dz 

"  Jro 

da 

*r« 

q’(Zra) 

Gri{i,p)  = 

-L  /  '  e« 

2ttj  ^rj 

G'ri(U)  = 

*  ' 

~t3l*dt 


(5.51) 

(5.52) 

(5.53) 

(5.54) 

(5.55) 


The  function  G<  is  expressible  in  terms  of  the  Incomplete  Airy  Functions 
given  in  Appendix  A.  To  apply  this  approximation  to  Jl5  the  proper  branch 
of  must  be  selected.  This  means  first  choosing  m  such  that 
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Figure  5.2:  Map  of  Re[g,  3(z)]  for  zn  -  z12  =  z,. 


0  ll/3  9  1/3 

— —  =  e-if+i'T;  m  =  0,±l,±2  (5.56) 

C(2>)J  <’(*■) 

agrees  with  the  arg  of  the  SDP  of  Figure  5.2.  The  indicated  choice  is  m  =  0 
so  that  arg (ht i)|,_0  must  equal  —  This  then  determines  the  choice  of 
branch  for  y/rji  since  it  must  be  consistent  with  Equation  5.52.  Then 


dz 

ds 
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so 


that 


c-iJ  =  e+J?{3-Jn) 

This  equation  is  satisfied  for  n  —  2  so  that 


l/a 

qi 

=  toil 

1/2  ( 

e+jtf 

qi 

=  toil 

e~: 

2. 

2-1, 1/1 

«11 

—  ft12 

g"(*n) 

1/2 


e  3  • . 


(5.58) 

(5.59) 

(5.60) 

(5.61) 


From  Equations  5.37,  5.38  and  5.48,  ara  ~  0  and  using  A. 28,  A. 31  and 
A. 37  Equation  5.46  becomes 

Ofcoo 


fco 0 


+  *  [/,(*„)!>„  -  01  (l<h I 0) 


fcl/3T7i 

fi(zu)hu  +/i(giz)hi2  _  /i(0) 

2qi  <?t(0)  J 


(5.62) 


Substituting  Equations  5.37  through  5.40  into  Equations  5.5  through  5.8 
and  making  use  of  Equation  5.61,  Equation  5.62  reduces  to 

/.(*)  ~  2ir/1(*11)/.„— e’»C; 


MiW 


+ 


Jb»/3 

[/i(2n)^n  _  /i(0) 
Vi  9i(°) 


(5.63) 


When  q  ^  0,  the  saddle  points  for  Ii  are  situated  as  Bhown  in  Figure  5.3. 

Turning  attention  now  to  /a,  the  proper  choice  of  arg(v/qj)  must  again 
be  made  and  m  must  be  selected  so  that 

1 1/3 

e-i f +/*?*.  m  =  0,  ±1,  ±2  (5.64) 


dz 

-2 

1/3 

2 

ds 

•=o 

c(».) 
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Figure  5.3:  Contours  of  Integration  and  Map  of  Re[g,(z)]  for  zn  ^  z13. 

agrees  with  the  arg  of  the  SDP  of  Figure  5.4.  The  indicated  choice  is  m  =  0 
so  that  arg(hji)|4=0  =  — The  choice  of  branch  for  y/r}[  is  then 

9  1/3  1/2 

n  =  0,  ±1,  ±2  (5.65) 

9  (*3i; 

so  that 

e-i*  =  e+if(1+2n).  (5.66) 
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This  equation  is  satisfied  for  n  —  —1  so  that 


Figure  5.4:  Contours  of  Integration  and  Map  of  Re[g,(z)| 
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/a(fc)  ~  ir[/2(zji)fe21  +  *G\  (-  |^|fcJ/3,o) 


+  [/a^ai^ai  —  /2(2aa)/ijal 

e*«(o) 


-fcao 


+ 


kllzr}\^7 

/a(2  21)^21  +  /a(2  22)^22  /a(0) 


(-  kal fca/3,0) 

(5.70) 


2»7a  ga(0)j 

which  using  Equations  5.37  through  5.40  and  Equations  5.5  through  5.8 
and  5.69  simplifies  to 

pkao 

h(k)  ~  2T/1(J„)i11-7jei!G;(-te|l’'3.0) 


pM°) 


lJfel/3 

/a(2ai)^ai  _  /a(0) 

9i(0)J 


*7a 


(5.71) 


Forming  the  sum  of  Ji  and  Ja  leads  to  further  simplifications.  Using 

9,(0)  =  9,(0)  (5.72) 

?!(0)  =  -9;(0)  (5.73) 

m  -  -t&  (5.74) 

(5.75) 

with  Equation  5.41  results  in 

m  =  /!(*)+/,(*) 

c*ao 

x  [c;  (+  M*"3,o)  +  g;  (-  |-h|  *3'3,o)]  (5.76) 


or 


I(k)  ~  2irf3(z,) 


x/a 


kV*q"(z.)\ 

X  [g;  (+  |>7.|*,/3,o)  +  g;  (-  m  *3/3,o)] .  (5.77) 
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While  Equation  5.77  is  complete,  it  is  not  in  the  best  form  for  comparison 
to  the  (singular)  first-order  stationary  phase  result.  By  rearranging  the 
solution  so  that  the  first-order  stationary  phase  solution  is  a  visible  factor, 
insight  may  be  gained  about  the  nature  of  the  corrected  solution.  The  usual 
first  order  isolated  saddle-point  solution  for  this  cylinder  is 

'<*>  ~  (s-78) 
Using  Equation  5.49,  Equation  5.77  may  be  rewritten  as 


m 


M**)t 


2ir 


*  |<(*#)| 


x  T(r}k7/3)  (5.79) 


where 


T(x)  =  [2v/irei*  |*|1/4e-J*l*|:’/a  (Ai*(+ |*|  ,0)  4- Ai*(- |x)  ,0))](5.80) 
f3  l3/3 

V  =  k(+2.) -9J(-«.)lj  (5-81) 

and  the  Airy  Functions  are  as  described  in  Appendix  A. 

From  the  above  equations,  it  is  seen  that  the  analysis  for  the  two  first- 
order  coalescing  saddle  points  results  in  a  form  that  can  be  written  as  the 
first-order  stationary  phase  evaluation  multiplied  by  a  function  T.  This 
function  acts  as  a  correction  to  the  singular  first-order  result.  Functions 
of  this  kind  in  UTD  are  called  Transition  Functions.  Since  it  is  known 
that  the  00  result  agrees  with  the  first-order  stationary  phase  evaluation 
for  the  reflected  field  far  away  from  the  pole,  it  is  anticipated  that  the 
transition  function  T(z)  is  the  desired  multiplicative  correction  to  GO.  The 
expected  behavior  of  the  transition  function  T(x)  then  for  large  argument 
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is  magnitude  of  unity  and  a  phase  of  zero.  This  is  seen  in  Figure  5.5. 
For  small  argument  z,  T(x)  should  approach  zero  in  such  a  way  that  the 
product  of  T(x)  and  the  infinite  GO  term  result  in  a  finite  and  accurate 
limiting  value  for  the  field  at  the  pole.  The  large  argument  form  of  the 
Incomplete  Airy  Function  is  given  in  Appendix  A.  Thus  for  z  0, 


For  e  «0, 


=  |*|1/4  «  1.453e*r/1>  |*|1/4 


(5.83) 


For  (0-f  9')  ~  0,  Eouations  4.1  and  4.5  show  that  the  order  of  the  singularity 
of  yfpc  is 

yfpc  «  [9  +  .  (5.84) 


For  v  =  3,  this  is  0(x  l^4).  The  behavior  of  ij  for  (0+9')  «  0  is  r\  oc  [0  4-  $']. 
Then  the  behavior  of  T(^Jfc,/,3)  is 


T(r,k7l*)<x\0  +  0')xlA .  (5.85) 

Since  this  is  a  zero  of  order  0(x+1^4),  then  the  GO  singularity  is  indeed 
cancelled. 
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Figure  5.5:  The  Transition  Function  T{z). 
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Phase  of  Transition  Function  (Degrees) 


5.5  The  Uniform  GO  (UGO)  Solution  for 

v  =  3 

The  total  corrected  GO  expression  for  the  reflected  field  from  a  superelliptic 

cylider  with  v  =  3  is  then 

e-ifc(p'+p) 

h;(p)  =  H0^fpc — — — 

where,  using  Equations  5.81,  5.10  and  5.36, 

2/3 

(5.87) 

and  </>,  0*,  pc  and  T(x)  are  given  by  Equations  4.1,  2.3,  4.6  and  5-80. 

For  small  (0  4-  0')/2,  rj  behaves  like 

(5.88) 

The  angular  region  where  the  pure  GO  result  is  invalid  may  be  determined 
via  Equation  5.87.  When  r/JkJ/3  >  10,  then  T(r]k3^3)  ~  1.  When  qfc*/3  <  3.5, 
then  T{r}k J/3)  begins  to  make  significant  corrections  to  the  GO  result.  A 
subjective  criterion  for  a  significant  departure  from  the  GO  result  (resulting 
in  a  greater  than  .04  dB  deviation)  is  rjk 2/a  <  5. 

The  numerical  results  UGO,  GO,  PO,  and  the  Method  of  Moments  are 
presented  in  Chapter  6. 

5.6  Reflected  Field  Behavior  for  v  ^  3 

The  correction  to  the  GO  solution  takes  the  form  of  a  two-parameter  mul¬ 
tiplicative  transition  function  T(x,  v)  where  v  is  the  superelliptic  “square- 


2  cos 


3/2\  1/3 


L  ( (« co®  (H£))3/3  +  (*“  (f¥i))3/J) 


X  T(r}k2/3)  (5.86) 
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ness”  parameter,  and  *  is  a  variable  describing  the  proximity  of  the  re¬ 
flection  point  to  the  pole,  or  zero-curvature  point.  The  transition  function 
T(x,i/)  does  two  things  near  the  pole;  it  acts  to  correct  the  GO  singularity, 
and  it  furnishes  an  aperture-like  oscillatory  pattern  behavior.  As  u  gets 
large  and  the  surface  around  the  pole  flattens,  this  aperture  effect  becomes 
more  and  more  pronounced.  In  the  limit  of  infinite  »/,  a  type  of 

pattern  is  expected. 

The  results  of  Chapter  5  explored  two  aspects  of  T(x,  v).  First,  the 
limiting  case  of  H ^  •  T(x  =  0,  u)  is  examined.  This  result  is  expressed 
in  terms  of  the  Gamma  Function  and  corresponds  to  an  evaluation  of  the 
field  reflected  from  the  pole  itself  for  arbitrary  u.  Second,  T(x,u  =  3)  is 
expressed  in  terms  of  Incomplete  Airy  Functions.  T(x,u  =  3)  is  associated 
with  the  reflected  field  both  near  and  far  from  the  pole,  but  only  for  the 
superellipse  whose  v  =  3. 

The  general  function  T(x,u)  cannot  be  expressed  in  terms  of  known 
functions,  but  it  is  hypothesized  that  the  general  transition  function  ac¬ 
counting  for  the  distributed  current  effects  around  the  pole  may  be  con¬ 
structed  heuristicaliy,  based  on  a  generalization  of  the  integral  form  of  the 
Incomplete  Airy  Function,  i.e. 

F(r}tu)  =  jf°  (e+in*  +  e“in‘)  e’$  dt  (5.89) 

=  2  /  coafafye**?  dt.  (5.90) 

Jo 

The  full  steepest-descent-path  analysis  is  not  tractable  for  the  general  case 
of  N  merging  saddle  points.  In  spite  of  this,  there  is  reason  to  believe  that 
Equation  5.90  may  form  the  basis  for  a  workable  T[x,  v)  after  all. 
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A  comparison  of  the  kernels  of  Equations  4.24,  5.10  and  5.90  reveal 
that  their  structures  are  the  same  in  the  neighborhood  of  t  =  0  for  all  u. 
It  is  not  therefore  unreasonable  to  suggest  that  Equation  5.90  is  in  fact 
the  correct  canonical  form  on  which  to  base  the  general  correction  T(x,  v). 
Further,  because  the  integration  increment  dt  is  strictly  a  real  quantity, 
Equation  5.90  easily  generalizes  to  all  real  v.  The  question  therefore  is 
whether  Equation  5.90  is  sufficiently  characteristic  of  superellipse  behavior 
to  construct  a  simple  and  accurate  T(x,*/),  or  whether  other  complicating 
factors  arise. 

As  v  becomes  large,  another  effect  appears,  which  is  quite  unrelated 
to  the  difficulties  associated  with  the  pole.  The  regions  where  the  radius 
of  curvature  decreases  with  increasing  v  develop  into  sharp  comers.  In 
these  regions,  GO  is  valid  provided  that  the  smallest  radius  of  curvature  is 
much  larger  than  the  wavelength.  A  natural  transition  from  the  mechanism 
of  reflection  over  to  diffraction  must  be  incorporated  into  the  solution  to 
achieve  full  generality  for  large  v  on  electrically  small  cylinders.  Some 
pertinent  results  concerning  reflection  by  surfaces  with  electrically  small 
radii  of  curvature  are  found  in  [4]. 
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Chapter  6 

NUMERICAL  RESULTS 


This  chapter  presents  the  numerical  data  for  the  reflected  field  from  vari¬ 
ous  superelliptic  cylinders,  calculated  by  the  various  methods  of  GO,  UGO, 
PO,  and  MoM.  First,  results  for  reflection  from  the  pole  are  presented,  then 
results  for  reflection  near  and  far  from  the  pole  for  v  —  3  cylinders  are 
presented.  The  curves  display  the  different  scattering  mechanisms  and/or 
artifacts  contained  in  the  various  methods.  Whenever  possible,  an  inter¬ 
pretation  of  the  results  will  be  pointed  out  in  the  discussions  accompanying 
the  graphs. 

With  regard  to  efficiency,  it  is  worth  noting  that  in  Table  6.1,  the  UGO 


a=b=lA 

a=b=3A 

a=b=5A 

a=b-10A 

a=b=30A 

GO 

.88  s 

.83  s 

.77  b 

.72  s 

.83  8 

UGO 

31.4  s 

16.9  s 

12.4  s 

8.25  s 

5.1  8 

PO 

36.1  s 

1.8  min. 

3.2  min. 

6.0  min. 

18  min. 

MoM 

1.7  min. 

4.8  min. 

9.5  min. 

32  min. 

>  4  hrs. 

Table  6.1:  Comparative  CPU  times  in  VPU  (VAX  780  Processing  Units) 
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solution  for  electrically  large  cylinders  is  more  than  two  orders  of  magnitude 
faster  than  Physical  Optics,  and  more  than  three  orders  of  magnitude  faster 
than  the  Method  of  Moments.  Further,  the  UGO  solution  is  the  only 
method  whose  computation  time  decreases  as  the  scatterer  gets  larger.  This 
phenomenon  is  due  to  the  asymptotic  behavior  of  the  transition  function 
T(x,  v  —  3);  it  is  easier  to  compute  for  large  arguments.  While  GO  is  also 
very  efficient,  it  does  not  produce  a  uniformly  valid  result  around  the  pole. 

The  slight  variations  visible  in  the  GO  CPU  times  reflect  the  variabilities 
and  inefficiencies  which  exist  in  a  multiuser  computer  environment.  The 
actual  computations  were  performed  on  a  VAX  8550  running  the  VMS 
operating  system.  Each  entry  in  Table  6.1  represents  90  backscattered  field 
computations  for  a  v  —  3  superquadric  cylinder. 

6.1  Reflected  Fields  from  the  Pole  for  Arbi¬ 
trary  v 

The  basic  asymptotic  result  concerning  the  backscattered  field  from  the 
pole  is  given  by  Equation  5.26.  This  equation  is  a  Physical  Optics  approx¬ 
imation  of  the  reflected  field  from  the  pole,  neglecting  the  second-order 
effect  of  creeping  waves,  and  excluding  false  PO  current  terminations.  The 
backscatter  field  for  various  types  of  superellipsoids  is  plotted  versus  v  in 
Figure  6.1.  Should  future  research  produce  a  general  reflection  transition 
function  T(x,t/),  Equation  5.26  will  be  useful  as  a  check  on  the  small  argu¬ 
ment  *  limit  of  T( 35,i/),  since  5.26  is  valid  for  all  v. 

Figure  6.1  illustrates  that  the  UGO  solution  provides  a  smooth  tran- 
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sition  from  (y^a/2)  type  of  behavior  which  is  characteristic  of  a  circular 
cylinder,  to  a  (2a)  behavior  which  is  representative  of  PO  scattering  from 
a  strip.  On  the  left  side  of  the  figure,  where  v  —  2,  the  backscattered  field 
of  a  circular  cylinder  as  a  function  of  radius  is  seen  by  the  intersection  of 
the  curves  with  the  leftmost  y-axis.  The  <Ja/2  behavior  is  evident.  On 
the  right  side  of  the  figure,  as  v  gets  large,  the  field  approaches  that  of 
the  backscattered  field  from  the  broad  side  of  a  rectangular  cylinder.  As  a 
function  of  radius,  the  intersection  of  the  curves  with  the  rightmost  y-axis 
shows  the  2a  behavior  expected  from  a  fiat  aperture.  Because  Equation 


5.26  is  only  valid  for  the  backscattered  field  from  the  pole,  Figure  6.1  does 
not  provide  information  about  the  field  pattern  away  from  the  pole. 
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BACKSCATTERED  FIELD  FROM  THE  POLE  (dB) 
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FAMILY  OF  CURVES 
OF  DIFFERENT  RADII 
A/A  -  B/A  - 
-1,2,3,4,5,6,7,8,9,10 


Figure  6.1:  Backseat tered  field  from  the  pole  for  different  radii  as  a  function 
of  v. 


Radius  -  b/X-a/X 


Figure  6.2:  Echo  width/7ra  from  the  pole  as  a  function  of  radius. 
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Normalized  Echo  Width  (dB) 


Figure  6.3:  Echo  width/jra  from  the  pole  as  a  function  of  radius. 
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In  Figures  6.2  and  6.3,  the  PO  results  of  Chapter  4  for  the  TE  polariza¬ 
tion  are  displayed  on  a  magnified  scale,  normalized  against  the  asymptotic 
results  of  Equation  5.26  and  Figure  6.1.  Since  the  data  were  generated 
using  numerical  integration,  the  endpoint  effects  of  the  false  PO  termina¬ 
tion  currents  make  themselves  visible  by  the  observed  periodic  oscillations. 
Note  that  the  endpoint  effects  diminish  both  as  the  radius  increases  and  as 
v  increases.  Also  note  that  the  scales  represent  0.2  dB  fluctuations  which 
accentuate  this  effect. 

The  decrease  of  the  endpoint  effects  with  increasing  radius  is  a  statement 
of  the  fact  that,  with  valid  assumptions  about  the  PO  currents  and  the 
surface  reflection  mechanism,  a  PO  result  approaches  GO  in  the  limit  of 
infinite  frequency.  Analytically,  PO  endpoint  contributions  usually  decrease 
as  0(l/fe),  where  the  stationary-phase  terms  are  typically  0(1). 

The  decrease  of  the  endpoint  effects  with  increasing  u  is  a  more  Bubtle 
effect,  but  is  easily  explained.  When  v  is  large,  the  cylinder  is  approxi¬ 
mately  a  rectangle  but  not  exactly.  The  current  terminations  do  not  occur 
at  the  corners,  as  in  the  cate  of  an  actual  rectangle,  but  the  constant-phase 
radiation  region  of  the  surface  does  end  there.  Instead  of  an  abrupt  current 
termination,  there  is  a  rapidly  fluctuating  phase  and  gradually  decreasing 
current  amplitude  on  the  top  and  bottom  sides  of  the  rectangle.  This  gen¬ 
erates  a  much  smaller  false  return  than  a  current  termination  immediately 
adjacent  to  the  constant-phase  large- amplitude  face  of  the  rectangle. 

The  validity  of  Equation  5.26  is  indirectly  confirmed  by  the  fact  that 
the  curves  for  all  v  tend  toward  zero  dB  with  increasing  radius.  If  Equation 
5.26  were  inaccuriie,  then  it  would  not  normalize  the  numerical  integration 
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to  unity  (0  dB).  Also,  if  all  endpoint  effects  were  subtracted  out,  then  all 
the  curves  should  lie  flat  on  the  zero  dB  line,  without  oscillation.  Again, 
it  must  be  noted  that  the  amplitude  scale  in  both  Figure  6.2  and  6.3  has 
been  greatly  exaggerated  to  display  the  termination  current  artifacts. 

0.2  Reflected  Fields  from  v  —  3  Superelliptic 
Cylinders 

This  section  presents  data  for  the  backscattered  fields  from  superquadric 
cylinders  of  v  =  3.  Here,  the  fields  are  investigated  for  reflection  points 
both  near  and  far  from  the  pole,  and  several  methods  are  plotted  together 
for  comparison.  Unless  otherwise  noted,  the  Method  of  Moments  results 
are  for  the  TM  polarization,  in  order  to  minimize  the  effect  of  creeping 
waves  around  the  cylinder.  Figure  6.5  is  the  exception,  which  shows  both 
the  TM  and  TE  results  using  the  Method  of  Moments. 

The  data  are  plotted  as  (echo  widths)/rra.  This  means  that  the  data  are 
normalized  to  the  2-D  echo  width  of  an  infinite  circular  cylinder  of  radius 
a.  Figure  6.4  shows  the  backscattered  field  as  a  function  of  the  angle  from 
the  x-axis  for  a  cylinder  of  radius  a/A  =  6/  A  =  1. 

The  GO  result  is  significantly  different  from  the  other  curves  in  virtually 
all  regions,  and  it  is  the  only  curve  that  is  unbounded.  Note  that  the  UGO 
and  the  PO  results  are  almost  indistinguishable  around  the  main  beam. 
The  UGO  and  the  Method  of  Moments  differ  by  less  than  1/2  dB  over  the 
entire  angular  range.  It  is  worth  noting  that  a/A  =  6/A  =  1  is  the  worst 
possible  case,  where  the  various  mechanisms  and  effects  such  as  creeping 
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I 
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waves  and  false  PO  returns  are  most  significant.  When  the  electrical  size 
of  the  cylinder  increases,  then  UGO  and  the  other  methods  converge.  This 
includes  GO  when  the  reflection  point  is  far  from  the  pole.  (Section  5.5 
and  Equation  5.87  define  “near”  and  “far”  from  the  pole.) 
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Figure  6.4:  Echo  width/wa  as  a  function  of  0  —  0',  (o/A  =  6/A  =  t). 
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Figure  6.5:  Echo  width/iro,  showing  both  TM  and  TE  polarizations. 
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Figure  6.6:  Echo  width/ira  a*  a  function  of  6  =  9',  (o/A  —  b/ A  =  3). 
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Normalized  Echo  Width  (dB) 


Figure  6.7:  Echo  width/rro  as  a  function  of  0  =  B\  (a/ A  =  b/\  —  5). 
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Normalized  Echo  Width  (dB) 


Figure  6.8:  Echo  width/jra  as  a  function  of  9  =  O',  (a/X  =  6/A  =  30). 
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In  Figure  6.5,  at  the  pole,  the  TM  polarization  via  Method  of  Moments 
is  1/2  dB  below  the  UGO,  whereas  the  TE  polarization  via  MoM  is  1/2  dB 
above  the  UGO  figure.  Slight  discrepancies  are  also  noticeable  around  the 
45  degree  mark.  The  TE  case  should  show  stronger  creeping  wave  effects 
for  small  radii.  This  could  explain  part  of  the  discrepancy.  Higher  order 
terms  in  the  asymptotic  analysis  could  also  account  for  these  effects.  In 
any  case,  the  comparisons  are  quite  close  considering  the  small  amount  of 
deviation,  thereby,  validating  the  results. 

The  differences  between  PO  and  UGO  in  Figures  6.4  and  6.6  -  6.8 
can  be  explained  by  the  termination  current  artifacts  present  in  the  PO 
solution,  which  are  absent  in  the  UGO.  The  UGO  solution  then,  is  the 
pure  reflection  mechanism,  inasmuch  as  the  PO  can  accurately  model  the 
primary  radiating  (stationary-phase)  currents  on  the  surface. 

As  the  radius  increases,  the  adverse  effects  diminish  and  the  UGO  result 
proves  to  be  very  accurate  indeed.  Figure  6.6  shows  that  for  a/ A  =  5/A  =  3, 
the  UGO  and  TM-MoM  results  differ  by  less  that  0.2  dB  over  the  entire 
angular  range.  Larger  radii  cases  in  Figures  6.7  and  6.8  show  even  closer 
agreement  between  the  various  methods. 
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Chapter  7 

SUMMARY  AND 
CONCLUSIONS 


Superquadric  surfaces  have  points  of  zero  curvature  at  which  GO  incor¬ 
rectly  predicts  an  infinite  reflected  field.  The  actual  reflected  field  is  finite 
and  is  well  approximated  by  the  method  of  Physical  Optics.  Rather  than 
simply  employ  PO,  however,  a  correction  to  the  GO  solution  is  constructed 
via  an  asymptotic  analysis  of  the  PO  formulation  for  the  reflected  field.  The 
purpose  of  the  asymptotic  analysis  is  to  achieve  an  improvement  over  PO 
by  avoiding  the  false  returns  associated  with  the  truncation  of  currents  at 
shadow  boundaries  and  by  avoiding  numerical  integration  over  electrically 
large  bodies.  In  the  spirit  of  UTD  then,  the  goal  is  to  retain  the  advan¬ 
tages  of  GO  in  terms  of  calculation  efficiency  and  analytic  simplicity  while 
simultaneously  enjoying  the  accuracy  and  physical  insight  afforded  by  the 
Physical  Optics,  hence  the  uniform  GO,  or  UGO  solution. 

The  numerical  results  of  Chapter  6  show  excellent  agreement  between 
the  UGO  solution  developed  here  and  the  Method  of  Moments.  There  is 
also  excellent  agreement  between  Physical  Optics  and  the  Method  of  Mo- 
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merits,  validating  the  inital  assumption  that  PO  is  a  good  way  to  describe 
the  reflection  mechanism.  In  most  of  the  regions  where  PO  and  the  Method 
of  Moments  disagree,  the  corrected  GO  solution  more  closely  agrees  with 
the  Method  of  Moments.  This  reveals  that  the  UGO  solution  avoids  the 
false  current  termination  effects  which  afflict  the  PO  result.  One  exception 
to  this  is  when  the  Method  of  Moments  includes  a  significant  higher  order 
scattering  mechanism  in  addition  to  the  reflected  field,  such  as  the  creeping 
wave.  In  this  case,  the  UGO  solution  does  not  include  this  effect. 

Finally,  the  issue  of  scattering  by  superellipses  in  three  dimensions  is 
suggested  as  an  interesting  area  for  future  research.  This  would  be  the  next 
natural  step  toward  making  superquadric  surfaces  into  a  viable  electromag¬ 
netic  modeling  tool. 
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Appendix  A 

COMPLETE  AND 
INCOMPLETE  AIRY 
FUNCTIONS 


The  Airy  Functions  satisfy  a  second-order  differentia)  equation  known  as 
Airy’s  Differential  Equation  A.2,A.3,A.23.  The  integral  solutions  of  these 
equations  are  known  as  Airy  integrals  and  their  properties  are  well  known. 
The  structure  of  the  Airy  integrals  is  such  that  they  are  characteristic  of  the 
arbitrary  interaction  of  two  saddle  points,  and  in  the  case  of  the  incomplete 
Airy  integral  they  describe  the  interactions  of  two  saddle  points  with  an 
endpoint  of  integration. 

This  appendix  presents  the  definitions  and  asymptotic  forms  of  the  com¬ 
plete  and  incomplete  Airy  functions.  While  the  Complete  Airy  Functions 
are  not  a  part  of  the  solutions  presented  in  Chapter  5,  they  are  used  in 
the  numerical  evaluation  of  the  Incomplete  Airy  Functions.  Asymptotic 
formulas  for  the  Airy  integral  are  given  in  (1). 


63 


A.l  Complete  Airy  Functions 

A. 1.1  Integral  Representation 


The  Complete  Airy  functions  Ai  and  Bi  are  defined  as 

Ai(i/)  =  ~  f  e^'-^dz  =  —  / 

2 nj  Jlx  2-k  Jli 

Bi(l)  =  si  .  =  4-  (  i. 

where  the  contours  Lnjz  are  shown  in  Figure  A.l. 


(A.l) 


Z-PIANE  S -PLANE 


Figure  A.l:  Contours  of  integration  for  the  Complete  Airy  functions 
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A. 1.2  Differential  Equations 


Ai "(z)  -  zAi(z)  =  0.  (A. 2) 

Bi "(z)  -  *Bi(2)  =  0.  (A.3) 

A. 1.3  Series  Representations 


Ai(z) 

=  Cif(z)  ~  c3g(z) 

(A.4) 

Bi(z) 

=  VZ[cif(z)  +  c2g(z)\ 

(A.5) 

/(*) 

“  to  V3A(3fc)! 

(A.6) 

,  1  3  1 '4  -  1 

=  1  +  3!Z  +~gTZ  +" 

•4-7  9 

9!  2  +- 

(A.7) 

M 

OO  /ox  -3M+1 

■  S'GW+u. 

(A.8) 

2  ,  2-5  7  2 

=r  ^  H - 2  3  4-  . .  2  4-  — 

4!  ^  7! 

•5;V+... 

10! 

(A. 9) 

“  +  i). 

=  i 

(A. 10) 

=  (3a  4-  l)(3a  +  4)  •  •  •  (3 

a  +  3Jb  -  1) 

(A. 11) 

Where  a  is  arbitrary  and  k  =  1, 2,3,  •  •  • 


d  =  Ai(0)  =  Bi(0)/>/3  =  3's/3/r(2/3) 

=  0.355028053887817  (A.12) 

c,  =  -Ai'(0)  =  Bi'(0)/Vl  =  3“,/3/r(l/3) 


0.258819403792807 


(A- 13) 


A.  1.4  Large  Argument  Forms 


Ai(*) 


Bi(*) 


for  | arg  z|  <  it, 


2-i/4  -  oo 


■  6  »* 


2y/n 
for  |arpz|  <  f, 
z-1/4 


fc=0 


E(-D‘ 


fc=o 


r(3fe  + 1) 

54fcfcir(fc  +  ±) 

r(3fc  +  l) 

54fcfcir(fc  + 1) 


(A. 14) 


(A. 15) 


A.1.5  Relations  between  Solutions 


Bi(z)  -  +  e'^Ai (ze~i3w/3)  =  0  (A.16) 

Ai(z)  +  ePKf*Ai(ze**'/9)  +  e~i,1r/3Ai(ze~i2ir/3)  =  0  (A.17) 

Bi(z)  +  ej2w/3Bi(zei3^3)  +  e~i7w/3Bi  (ze~i7^3)  =  0  (A.18) 

2  Ai(ze±ia,r/3)  =  e±ia*/3  [Ai(z)  ±  Bi(z)]  (A.19) 


66 


A.2  Incomplete  Airy  Functions 


The  Incomplete  Airy  functions  are  defined  as 

=2^  e^*+V)ds,  i  =  1,2,3 

(A. 20) 

1  fOOjei+i  j 

=  —  /  e^-^dz 
2  irj  Jjfii 

0  <  V»i  <  f,  -*<*<-*  (A. 21) 

where  the  ooe-^  and  /?,■  correspond  to  the  endpoints  of  the  S-plane  paths 
Lj  shown  in  Figure  A.2.  The  contour  Li  may  be  deformed  onto  the  positive 


by  Ai 

Ai(,,0)  s  0.(*«  =  l  £  «*-*>*  =  JL  j‘~  (A.22) 

A. 2.1  Differential  Equation 


^  .  a ' 

a*a  2  Ja/? 


XI(z,/3)  =  0 


A.2.2  Large  Argument  Forms 


(A. 23) 


Assume  that  /9  >  0  and  *  >>  0. 

j  eX\P3+*0) 


Ti(z,0) 

■SI(-x,/3) 

where 


2n  (32  +  x  ’ 

A.-jj*"* 

2»r 


/?*  +  x  >>  0 


-g(°) 


(A. 24) 


(A. 25) 


5(«) 

✓(«) 

</"(*) 

P(0) 

y(o) 

5"(0) 


dt  __  2s 
da  t2  —  * 

7  [<K*)  -  <03(*)] 

“j  [3tV(«)  -  3<53(j)  -  sy4(s)] 

J_ 

x4 

J_ 

_  3x 
5 

12e7/4 


(A. 26) 
(A.27) 


*  -  (^4’f 

a  —  ±sa  for  db  (/?  -  x1/2) 
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A. 2.3  Relations  between  Solutions 


G,(v,0)  =  G,(V,/3)  -  Ai(V)  (A.28) 

G3(v,0)  =  G1(^/?)-|[Ai(t?)+jBi(,7)]  (A. 29) 

=  GM)  +  \[ Ai(r?)  -  jBi(^»  (A.30) 

Ai(ij)  =  G1(r)if3)  +  Gl(v\0)  (A.31) 

Gi(Tje+’2*f3,0)  =  e~*2ir/3G2(tu0e+i2w'3)  (A.32) 

Gi(r}e~*,"t3,(3)  =  e+**'3G3{v,0e-i2w/3)  (A.33) 

G3(ve+ill,,3,0)  =  e~»*'3G3{i h(3e+i2*'3)  (A.34) 

G2(Ve~i2*/3t0)  =  e+^Gt^fle-i2*'3)  (A.35) 

G3(rfe+i2n/3,0)  =  e-i2*/3Gi{ri,l3e+i2w'3)  (A.36) 

G3{Ve~^3,0)  =  e+iu'3G2{T,}0e'i2''3)  (A.37) 
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Appendix  B 

COMPUTER  PROGRAMS 


This  appendix  contains  the  various  computer  programs  used  to  calculate 
the  reflected  fields  from  the  superelliptic  cylinder.  Some  of  the  Airy  func¬ 
tion  subroutines  were  contributed  by  M.  C.  Liang. 


PROGRAM  UNIFORM  GO 

This  program  calculates  the  reflected  field  from  the  superellipse  using  Uni¬ 
form  GO,  or  UGO. 

options  /extsnd-source 
program  go 
I 

!  This  program  calculates  the  00  backscattered  field  from  a 

!  superellipse . 

! 

parameter  pi  *  3.141592664 

real  theta,  r_c,  Hgo,  a,  b,  knu,  area,  Refpnt,  zs,  zc,  t,  tp 


real  x,  xxx 
complex  transition 


external  r_c,  transition 


Ref pnt (theta)  *  atan(sign(abs( 

+  ((b/a)**knu)*tan(theta))**(l/(knu-l)) ,tan(theta))) 

Hgo  (  theta  )  ■  sqrt  (  r_c(  Ref  pnt  (theta) ,  a,  b,  knu  )/2  ) 

knn  *  3.0 

type*,  ’Input  a’ 
accept* ,  a 
b*a 

do  theta*-46,  46,  .995 

Calculate  [  x*sigma*k?2/3)  ] 

t  *  aba (theta  *  pi/180) 
tp  =  t 

zs  =  (  b  *  sin((t+tp)/2)  )  **  (3./2.) 
zc  *  (  a  *  cos((t+tp)/2)  )  **  (3./2.) 

xxx  *  Z8**(2./3.)  *  (  2*cos((t-tp)/2)  / 

+  (  zs  +  zc  )**(l/3.)  )**(2/3.) 

xxx  *  ain  (  xxx*(2*pi)**(2/3.) ,  24.0  )  !  bug  in  airy  function 

!  for  large  arguments 

area  ■  Hgo(  theta*pi/180  )  !  *  cabs  (  transition(x)  ) 
area  *  area* *2  *  2*pi  /  (  pi  *  a  ) 
write(6,*)  theta,  area 
end  do 

end 
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SUBROUTINE  RC 

This  subroutine  calculates  the  radius  of  curvature  at  a  point  on  the  surface 
of  a  superelliptic  cylinder. 


options  /extend-source 

real  function  R-C(  phi,  a,  b,  knu  ) 

real  phi,  s,  c,  knu,  a,  b 

real  sa,  ca,  snail,  caml,  sam2,  cam2 


s  3  abs(  sin(phi)  ) 

c  =  abs(  cos (phi)  ) 

sa  sB*4i  knu 
ca  *  c  *•  knu 

saml  =■  s  **  (knu-1) 

caml  *  c  **  (knu-1) 

sam2  *  8  ++  (knu- 2) 

cair2  «  c  *+  (knu-2) 

R_C  *  (  ((a*+knu)*saml)**2  +  ((b**knu)*caml)**2  )**(3./2.)  /  ( 


+  (knu-1)  * 

+  (a*b)**(knu-l)  *  (  (a*s)+*knu+(b*c)**knu  )**(l+l/knu)  * 

+  (sec)  **(knu-2)  ) 

return 

end 


FUNCTION  TRANSITION 

This  subroutine  calculates  the  value  of  the  transition  function  T(x). 
options  /extend-source 
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complex  function  Transition  (  argument  ) 

complex  iai,  j,  z 

external  iai 

real  pi ,  argument 

parameter  (  pi  ■  3.141592654,  j»(0.,l.)  ) 

*! 

z  «  iaiC  0.0,  + argument  ) 

+  +  iai(  0.0,  -argument  ) 

2  »  conjg(z)  *  2  *  aqrt(pi)  *  (abs (argument )**(1 ./4. )) 

+  *  exp(  -j*2./3.*(aba(argument)**(3./2.))  ♦  j*pi/4  ) 


Transition  ■  z 

return 

end 
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PROGRAM  POLE  BACKSCATTER 

This  program  calculates  the  backscattered  field  from  the  pole  of  a  superel- 
liptic  cylider. 


program  super  ellipse  pole 
integer  n 

real  k,  c,  a,  x,  y,  pi,  gamma,  btan2,  knu,  fact,  deriv 
complex  i,  j/(0.,l.)/ 

pi  =  3.141692654 
k  «*  2  *  pi 

type*,  ’input  a’ 
accept*,  a 

do  knu»2.0,  20.0,  .1 

deriv  ■  -  2  *  gamma (  knu  )  i  -2(alpha-l)! 

i  ■  gamma(l/knu)  *  2/knu  * 

♦  (-gamma(knu+l)/(deriv*j*k*a))**(l/knu) 

+  *  cexp(j*k*2*a) 

+  *  C8qrt(j*k/(2*pi))  *  a 

WRITE  (6,*)  knu,  cabs(i),  btan2 ( aimag ( i ) ,real(i))*180/pi,  i 


end  do 
END 


real  function  gamma(x) 

integer  i 

real  x,  fact 
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call  gmiiwm  (  x,  fact, 
gamma  ■  fact 
return 
end 


PROGRAM  PHYSICAL  OPTICS 

This  program  was  used  to  calculate  the  Physical  Optics  reflected  field  from 
the  superelliptic  cylinder. 

OPTIONS  /EITEND.SOURCE 

PROGRAM  BISTATIC  PO 

INTEGER  PIECES,  I,  N,  MAIPNTS 

real  btan2 

external  btan2 

C0MPLEI*16  C,  J,  H,  z 

real  phase 

REAL*8  PHI,  LOW,  HIGH,  RE_H,  IM-H,  RE-RESULT,  IM-RESULT,  lambda 
REAL*8  FMIN,  FMAI,  FSTEP ,  F 

REAL*8  PI,  A,  V,  K,  L,  U,  H.N0RM,  i.-.in,  amax,  NORM-H ,  ONE,  CDABS 

COMMON  A,  PHI,  J,  K,  V 
EXTERNAL  RE-H,  IM-H 

phase(z)  ■  btan2(  dimag (2),  real(dreal(z))  )*180/pi 

ONE  -  1.0000000000000000 
J  *  (0. ,1.) 

PI  *  4  *  ATAN  (ONE) 

TYPE*,  ’Input  the  observation  angle  (degrees)' 

ACCEPT*,  PHI 

TYPE*,  ’Input  the  cylinder  radius  (meters)’ 

ACCEPT*,  A 

TYPE*,  ’Input  the  frequency  range  [low, high .step] (Ghz) ’ 

ACCEPT*,  fain,  fmax,  fstep 

TYPE*,  ’Input  the  superelliptic  parameter  (knu) ' 

ACCEPT*,  v 


phi  *  phi  *  pi/180 
fmin  ■  fain  *  le9 
fmax  •  foaz  *  le9 
fstep  *  fstep  *  le9 

i 

!  Loop  through  several  frequencies. 

j 

maxpnta  *  nint(  (fmar-fmin) /fstep  ) 

WRITE(3 ,*)  maxpnta,  real (fmin/ la6) ,  real (f atop/ le6) 
DO  N*l,  M11PNTS  +  1 

; 

!  Change  aoma  parameters  for  each  point . 

i 

C  A  -  (  (M- 1 ) * ( amax-aain) /MAXPNTS  ♦  amin  ) 

C  PHI  «  (  (N-l)*(pi/2)/MAXPNTS  +  -pi/2  ) 

j 

f  -  (  (*f~  1 )  *  (f aax-foin)  /MA1PNTS  +  fain  ) 
lambda  *  299  792  837.1  /  f 
k  ■  2  *  pi  /(  lambda  ) 
c  ■  2  *  CDSQRT((j*k)/(8*PI)) 

j 

!  Calculate  the  point. 

; 

HIGH  -  PI/2+ATANC  ABS(TAN(PHI))**(1/(V-1)) 

+  *SIGN(-ONE,SIN(PHI))  ) 

LOW  -  HIGH  -  PI 
H  ■  (0.,0.) 

PIECES  -  2  *  NINT(  A/lambda  ) 

DO  1-1,  PIECES 

i 

L  «  (1-1)  *  ((HIGH  -  LOW) /PIECES)  +  LOW 
U  -  I  *  ((HIGH  -  LOW) /PIECES)  +  LOW 

t 

!  Integrate  the  real  and  imaginary  part a  of  H  due  to  J. 

! 

CALL  DQG32(  L,  U,  RE_H,  RE-RESULT  ) 


CALL  DQG32(  L,  U,  IM_H,  IM.RESULT  ) 
H  «  H  ♦  CMPLKRILRESULT ,  IMJRESULT) 


i 

END  DO 

H  *  H  *  C  *  A 

!  Normalize  the  result  against  the  GO  result. 

N0RM-H  -  A  *  PI 

ELNORM  »  2  *  PI  *  CDABS (  H  )**2  /  NORM-H 

!  Output  the  data  to  a  file  and  notify  the  user  of  progress. 

WRITE(3,*)  20*logl0(cdabs(H)),  phase(h) 

TYPE*,  N,  »  of  *,  MAIPNTS,  *  f«’,  f/le9,  ’  h«’, 

+  20*logl0(real(abs(h_norm))) 

i 

END  DO 

CALL  STATISTICS 
END 


These  functions  compute  the  real  and  imaginary  parts  of  H. 


! 

f 


REAL*8  FUNCTION  RE-H  (  I  ) 

REAL*8  I,  AMP-NUMER,  PHASE.NUMER,  PHASE-DENOM,  AMP_DEN0M 

REAL*8  PHI,  A,  S,  C,  V,  TMP,  IM_H,  K 

COMPLEX* 16  J 

COMMON  A,  PHI,  J,  K,  V 

LOGICAL  REAL  /.FALSE./ 

REAL  *  .TRUE. 


!  The  real  and  imaginary  parts  are  very  similar. 


ENTRY  IM-H  (  X  ) 


i 
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S  =  SIN(X) 

C  »  COS(X) 

Calculate  the  numerators  of  the  amplitude  and  phase  functions. 

AMPJfUMER  «  SIN(PHI)*(SIGN(ABS(S)**(V-1),S))  + 

+  COS (PHI) ♦(SIGN(ABS(C) **(¥-!) ,C)) 

PHASE JfUMER  ■  2  *  (  C  ♦  COS (PHI)  +  S  *  SIN (PHI)  ) 

Build  the  denominators  of  the  amplitude  and  phase  functions. 

TMP  *  (  ABS(S)**V  +  ABS(C)*+V  ) 

PHASE-DENOM  «  TMP**(1/V) 

AMP-DENOM  -  PHASE-DENOM  *  TMP 

Now  combine  the  functions. 

IF  (  BEAL  )  THEN 
RE-H  -  ( AMPJfUMER/ AMP-DENOM)  * 

+  C0S(  K*A*  (PHASE  JIUMER/PHASE-DENOM)  ) 

ELSE 

IM_H  •  (AMP JTUMER/ AMP-DENOM)  * 

+  SIN(  K* A* (PHASE-NUMER/PHASE-DENOM)  ) 

END  IF 

Asume  imaginary  unless  the  real  entry  point  was  taken. 

REAL  -  .FALSE. 

END 
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SUBROUTINE  CDQG32 

This  subroutine  performs  a  complex  double-precision  Gaussian  quadrature 
numerical  integration  on  a  given  function. 

+  !  + 

* ! 

* !  This  is  a  complex  double-precision  32-point  gaussian  quadrature 
*!  integration  routine.  It  integrates  the  function  FUNC  on  a  strait 

*!  line  in  the  complex  plane  from  ZL  to  ZU. 

*! 

*!- 

DOUBLE  COMPLEX  FUNCTION  CDQG32  (  ZL,  ZU,  FUNC  ) 

DOUBLE  COMPLEI  FUNC,  FCT,  Y,  ZL,  ZU 

DOUBLE  PRECISION  A,  Cl,  C2,  C3,  C4,  C5,  C6,  C7,  C8,  C9,  CIO 

DOUBLE  PRECISION  Cll,  C12,  C13,  C14,  C15,  C16,  C17 
DOUBLE  PRECISION  X,  Bl,  B2,  B3,  B4,  B6,  B6,  B7,  B8,  B9,  BIO 


DOUBLE  PRECISION 
PARAMETER  (  A 
PARAMETER  (  Cl 
PARAMETER  (  C2 
PARAMETER  (  C3 
PARAMETER  (  C4 
PARAMETER  (  C6 
PARAMETER  (  C6 
PARAMETER  (  C7 
PARAMETER  (  C8 
PARAMETER  (  C9 
PARAMETER  (  CIO 
PARAMETER  (  Cll 
PARAMETER  (  Cl 2 
PARAMETER  (  C13 
PARAMETER  (  C14 
PARAMETER  (  C1B 


BU,  B12 ,  B13,  B14,  B16,  B16,  B17 
0.5D0  ) 

.49863193092474  DO  ) 
.49280576677263  DO  ) 
.48238112779375  DO  ) 
.46746303796886  DO  ) 
.44816067788302  DO  ) 
.42468380686628  DO  ) 
.39724189798397  DO  ) 
.36609105937014  DO  ) 
.33152213346510  DO  ) 
.29385787862038  DO  ) 
.25344995446611  DO  ) 
.21067563806631  DO  ) 
.16693430114106  DO  ) 
.11964368112606  DO  ) 
.07223698079139  DO  ) 
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* ! 

*? 


PARAMETER 

( 

C16 

m 

.024153832843869  DO 

) 

PARAMETER 

( 

B1 

m 

.35093060047350 

D-2 

) 

PARAMETER  ( 

B2 

m 

.8137197386462 

D-2  ) 

PARAMETER 

( 

B3 

m 

. 12696032654631 

D-l 

) 

PARAMETER 

< 

B4 

a 

.17136931456510 

D-l 

) 

PARAMETER 

( 

B5 

m 

.21417949011113 

D-l 

) 

PARAMETER 

( 

B6 

a 

.25499029631188 

D-l 

) 

PARAMETER 

( 

B7 

8 

.29342046739267 

D-l 

) 

PARAMETER 

( 

B8 

8 

.32911111388180 

D-l 

) 

PARAMETER 

( 

B9 

8 

.36172897054424 

D-l 

) 

PARAMETER 

( 

BIO 

8 

.39096947893535 

D-l 

) 

PARAMETER 

( 

Bli 

8 

.41666962113473 

D-l 

) 

PARAMETER 

( 

B12 

8 

.43826046502201 

D-l 

) 

PARAMETER 

e 

B13 

8 

.45586939347881 

D-l 

) 

PARAMETER 

( 

B14 

8 

.46922199640402 

D-l 

) 

PARAMETER 

( 

B15 

8 

.47819360039637 

D-l 

) 

PARAMETER 

( 

B16 

a 

.48270044267363 

D-l 

) 

FCT  (  X  ) 

a 

FUNC 

ZL  +  (ZU-ZL)*X 

) 

Y 

8 

( 

0 

0.  ( 

).C 

) 

Y 

8 

Y 

+ 

B1 

4 

( 

FCTU+Cl) 

4 

FCT(A-Cl) 

) 

Y 

8 

Y 

+ 

B2 

* 

( 

FCTCA4C2) 

4 

FCTU-C2) 

) 

Y 

8 

Y 

+ 

B3 

4 

( 

FCTU+C3) 

4 

FCTU-C3) 

) 

Y 

8 

Y 

+ 

B4 

* 

( 

FCT(A4C4) 

4 

FCTU-C4) 

) 

Y 

8 

Y 

+ 

B5 

4 

( 

FCTU+C6) 

4 

FCTU-C5) 

.) 

Y 

8 

Y 

+ 

B6 

* 

( 

FCT(A4C6) 

4 

FCTU-C6) 

) 

Y 

8 

Y 

+ 

B7 

4 

( 

FCT(A+C7) 

4 

FCTU-C7) 

) 

Y 

8 

Y 

+ 

B8 

* 

( 

FCTCA+C8) 

4 

FCTU-C8) 

) 

Y 

8 

Y 

4 

B9 

* 

( 

FCT(A4C9) 

4 

FCTU-C9) 

) 

Y 

8 

Y 

4 

BIO 

* 

c 

FCT(A+C10) 

4 

FCT(A-CIO) 

) 

Y 

8 

Y 

+ 

Bll 

« 

( 

FCT(A4Cli) 

4 

FCT(A-Cll) 

) 

Y 

8 

Y 

+ 

B12 

4 

( 

FCTU+C12) 

4 

FCTU-C12) 

) 

Y 

8 

Y 

+ 

B13 

4 

( 

FCT(A4C13) 

4 

FCTU-C13) 

) 

Y 

8 

Y 

+ 

B14 

4 

( 

FCT(A+C14) 

4 

FCTU-C14) 

) 

Y 

8 

Y 

4 

B16 

4 

( 

FCTU4C16) 

4 

FCT(A-C16) 

) 

Y 

a 

Y 

4 

B16 

4 

( 

FCTU+C16) 

4 

FCTU-C16) 

) 
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CDQG32  «  (ZU-ZL)  *  Y 

RETURN 

END 
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FUNCTION  INCOMPLETE  AIRY 

These  functions  compute  values  for  the  Incomplete  Airy  functions. 


complex  function  ini  (  b,  x  ) 
complex  iaip,  ibi,  ibip,  zl,  z2,  z3,  z4 
real  b,  x 

call  aiinc  (  b,  x,  zl,  z2,  z3,  z4  ) 

iai  ■  zl 

return 

entry  iaip 

call  aiinc  (  b,  x,  zl,  z2,  z3,  z4  ) 

iaip  «  z2 

return 

entry  ibi 

call  aiinc  (  b,  x,  zl,  z2,  z3,  z4  ) 

ibi  »  z3 

return 

entry  ibip 

call  aiinc  (  b,  x,  zl,  z2,  z3,  z4  ) 

ibip  ■  z4 

return 

end 

SUBROUTINE  AIINC (BETA , IS , All .AIIP.AIH, A1HP ) 

COMPLEX  All , AIIP , FCT , FIST , C J , CJP4 
COMPLEX  AIH , AIHP , ATEM (1001) .ATMP(lOOl) 

COMPLEX  61 , GIP , HI , HIP , AI , AIP , BI , BIP 
COMPLEX  CT1 , CT2 , CT3 , CT4 , CT5 
C0MM0N/ERR/ERR1 , ERR2 , ERR3 , ERR4 
C 

C  APRIL  17,  1988 

C  ************•****•*************+******************************* 
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C  *  THIS  SUBROUTINE  IS  USED  TO  CALCULATE  THE  INCOMPLETE  AIRY  * 

C  *  FUNCTION  AND  ITS  DERIVATIVES,  THE  FUNCTION  IS  DEFINED  AS  * 


C  *  AII(B,S)«INT(B,INF,FS)  * 
C  *  AIIP (BS) *INT(B , INF , C J+T*FS )  * 
C  *  FS»EIP(CJ*(T**3/3+S*T))/TPI  * 
C  *  NOTE  THAT  THE  AIRY  FUNCTION  AI  IS  DEFINED  AS  * 
C  *  AI(S)=INT(“INF, INF.FS)  * 


C  ********** ******************* **************************** ****** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IARG-l :  SMALL  NEGATIVE  ARGUMENT  OR  POSITIVE  ARGUMENT  FORM; 
WITH  BETA  »1 

-2:  LARGE  NEGATIVE  ARGUMENT  FORM(  ZS  <<0); 

USING  THE  FRESNEL  INTEGRAL  (ASYMPTOTIC  FORM) 

*3:  SMALL  NEGATIVE  ARGUMENT  OR  POSITIVE  ARGUMENT  FORM; 
WITH  0<  BETA  «1 

-4:  EVALUATE  THE  INTEGRAL  FORM;  AIINC(B,S)»AIINC(0,S) 
-INT(O.B.FS);  WHERE  FS*EIP(CJ*(0.33*T*T*T+S*T))/TPI 


ICQM-1:  IF  BETA  <  0  ;  IN  THIS  CASE  TAKE  THE  COMPLIMENTARY 
PART  OF  AIKBETA.XS); 

I.E.  AII-AI(XS)-AII_*_(ABS(BETA)  ,IS) 

■0:  IF  BETA  >  0  . 

PI-3. 14169266369 
TPI-2 . *PI 
CJ-(0.,1.) 

CT1— CJ*PI/4. 
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CJP4-CEIP (CT1) 

QPI-SQRT(PI) 

BETl-ABS(BETA) 

ICOM-O 

IF (BETA . LT . 0 . ) ICOM* 1 
IARG-4 

T1-BET1*BET1+IS 
IF(T1.GE.18.)IARG»1 
C  IF (BET1 . LE , ERR4) IARG*3 

IF (IS . LE . -7 . ) IARG*2 

CALL  GIHI (IS , GI , GIP , HI , HIP , AI , AIP ,BI , BIP) 

AIH-0.5*(AI+CJ*GI) 

AIHP-0 . 6* (AIP+CJ*GIP) 

C 

GO  TO  (100 ,200 ,400, 400) IARG 
100  CONTINUE 

C 

C  THIS  IS  THE  LARGE  BETA  FORM;  WITH  SMALL  NEGATIVE  IS  OR 

C  LARGE  POSITIVE  IS 

C 

CT1-CJ* (BETl+BETi*BETl/3 . +IS*BET1) 

CT2*CEIP(CT1) 

BBS 1*1 . /(BET1*BET1+IS) 

BBS2-BBS1*BBS1 

BBS4«BBS2*BBS2 

CT3»CJ*BBS1* ( 1 . - ( 10 . *BETl*BETl-2 . *IS) *BBS4) 

C  CT4*BBS1*BBS2*2 . *BET1* ( 1 . -16 . * (3 . *8ET1*BET1-IS) *BBS4) 

CT4*BBS1*BBS2*2 . *BET1 
All- ( CT3+CT4) *CT2/TPI 
C 

C  CT3-CJ*BBS2* (-1 . +BBS4* (62 . *BETl*BETl-8 . *IS) ) 

C  CT4-BBS4*8 . *BET1*( - 1 . + ( 1 10 . *BET1*BET1-30*IS) *BBS4) 

C 

CT3—CJ*BBS2 
CT4— BBS4*8 .  *BET1 

AIIP* (CT3+CT4) *CT2/TPI+CJ+BET1 *AII 
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GO  TO  900 
200  CONTINUE 

C 

C  THIS  IS  THE  ASYMPTOTIC  FORM  FOR  LARGE  NEGATIVE  IS; 

C  IN  THIS  CASE  THE  FUNCTION  IS  APPROXIMATED  BY  FRESNEL 

C  INTEGRAL 

C 

C  THE  LARGE  ARGUMENT  IS  DECIDED  TO  BE  USED  FOR  XS<  -7 . ; 

C  IN  THIS  CASE  THE  ERROR  IS  ABOUT  0.3X  FOR  BETA»0. 

C 

ABXS-ABS(XS) 

SAIS-SQRT(ABIS) 

SSAI-SQRT(SAIS) 

DEPS-BET1-SAXS 
CT3—C  J*2 .  *ABXS*SAXS/3 . 

CT4«CEXP(CT3) 

C 

C  BRANCH  INTO  LARGE  AND  SMALL  ARGUMENT  FORMS 

C 

IF(ABS(DEPS) .GT.O .O3)G0  TO  250 
C 

C  SMALL  ARGUMENT  FORMS 

C 

Tl-1 . +DH>S/ (6 . *SAXS) -DEPS*DEPS/ (72 . *ABIS) 

DEBTA-T1*DEPS*SSAI 

DEBTA2«DEBTA*DEBTA 

CT1-CJ+DEBTA2 

CT2-CEXP(CT1) 

CT1«QPI*(0. 7070108781, 0.7070106781) 

FINT*0 . 6*CTl-DEBTA-0 . 33333334*DEBTA*DEBTA2*CJ 
Tl»l . -0 . 625*DEPS/SAIS+0 . 3402777778*DEPS*DEPS/ ABIS 
CT5—0 . 6*CT2*T1*CJ/  (ABXS*3 . ) 

AII-CT4* ( (1 . +CJ*0 . 10418666867/ (SAXS+ABXS) ) *FINT/SSAI+CT6) /TPI 

T2"l . -0.4375*DEPS/SAXS+0.409722222*DEPS*DEPS/ABXS 
CT6-T2*CT2/(3.*SAIS) 

AIIP-CT4* ( <C J+0 . 1458333333/ ( SAXS*ABIS) ) *SSAX*FINT~CT6) /TPI 
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GO  TO  900 


C 

C  LARGE  ARGUMENT  FORMS 

C 

260  GEO-l./SSAX 

Tl=BETl+BETl*BETl/3 . +XS*BETl+2 . *SAIS* ABIS/3 . 

DEBTA-SQRTCT1) 

IF  (BET1 .  LT .  SAXS  )  DEBTA— DEBTA 
C  GEBTA-GEO 

C  T1»BET1*BET1+XS 

C  IF ( ABS (T1 ) . GE . 0001) GEBTA=2 . *DEBTA/T1 

GEBTA-2 . *DEBTA/ (BET1*BET1+XS) 

DEBTA2=DEBTA*DEBTA 
CALL  FCTX ( 1 , FCT , DEBTA2) 

CT1=CJ*DEBTA2 

CT2*CEXP(CT1) 

CT1-QPI* (0 . 7070106781 , 0 . 7070106781) 

C 

FINT-0 . 5*C0NJG(FCT) *CJ*CT2/DEBTA 
IF (DEPS . LT . 0 . )FINT«CT1+FINT 
C 

CT6-0 . 5*CJ* (GEBTA-GEO) /DEBTA*CT2 

AII-CT4* ( (GEO+C J  *0 . 10416666667*SSAX/ ( ABXS* ABXS) ) *FINT+CT5) /TPI 
C 

CT6- (BET1/ (BET1*BET1+XS) -0 . 5+SSAX/DEBTA) *CT2 

AIIP-CT4* ( (C J*SSAI+0 . 1458333333/ (SSAX*ABIS) ) *FINT-CT5) /TPI 

GO  TO  900 
C 

C300  CONTINUE 
C  BET2-BETA+BETA 

C  BET4-BET2*BET2 

C  BET8-BET4*BET4 

C  XS2-IS*XS 

C  XS4-IS2*XS2 
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C  T1=BET1-BET1*BET2* (BET4/126 . +XS*BET2/15 . +XS2/6 . ) 

C  T2*BET2* (BET2/ 12 . +0 . B*XS) - (BET4/ 12 . ) * (BET4*BET2/136+ 

C  $  XS*BET4/12 . +XS2*BET2/3 . +IS+XS2/2 . ) 

C  T3= (BET1*BET4/ 6 . ) * (BET8/4212 . +IS*BET2*BET4/297 . 

C  $  +IS2*BET4/S4 . +XS*XS2*BET2/21 . +XS4/20 . ) +T1 

C  AII*AIH-CMPLX(T3,T2)/TPI 

C  T4»- (BET1+BET2/3 . ) * (BET2/5 . +XS) ♦ (BETl*BET4/6 . ) 

C  $  *(IS*XS2/6.+XS2*BET2/7.+XS*BET4/27.+BET2*BET4/297.) 

C  T6*0 . 6*BET2* ( 1 . -BET2*BET4/72 . -XS*BET4/9 . -XS2*BET2*0 . 26) 

C  AIIP*AIHP-CMPLX(T4,T5) /TPI 

C  GO  TO  900 

400  CONTINUE 

N3«26 

IF((BET1.GE.3. ) .AND. (XS.GE.O. ))N3=60 

NNT=BET1+N3+1 

Nl*NNT/2 

N2=«NNT-N1*2 

IF ( N2 . EQ . 0) NNT*NNT+ 1 

IF (NNT . LE . 3)NNT*3 

DELT«=BET1  /  (NNT- 1 ) 

C 

DO  410  I-l.NNT 
T1«DELT*(I-1) 

CTl®CJ*(Tl+Tl*Tl/3.+XS+Tl) 

CT2=CEXP(CT1) 

CT3«CJ+T1*CT2 
ATEM(I)«CT2 
ATHP(I)-CT3 
410  CONTINUE 
C 

Nl**(NNT-i)/2 
CT4«(0. ,0.) 

CT6-(0. ,0.) 

DO  440  J«1,N1 
J1»2*J 
J1M-J1-1 
J1P-J1+1 
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CT4-CT4+ ( ATEM ( JIM) +4 . * ITEM ( J 1 ) +ATEM ( J IP) ) *DELT/3 . 

CT5-CT5+ (ATMP (JIM) +4 . *ATMP ( J1 ) +ATMP ( J1P) ) *DELT/3 . 

440  CONTINUE 
C 

AII=AIH-CT4/TPI 
AIIP-AIHP-CT6/TPI 
900  IF (ICOM.EQ.C) RETURN 

C 

C  ICOM-O  MEANS  THAT  BETA  >  0;IC0M*1  MEANS  THAT  BETA  <  0. 

C 

AII=AI-CONJG(AII) 

AIIP«AIP-CONJG(AIIP) 

RETURN 

END 

SUBROUTINE  GIHI (IS , GI , GIP , HI , HIP , All , AI1P , BI1 , BI1P> 

C 

C  XS:  REAL  ARGUMENT 

C  GI, GIP .HI, HIP:  COMPLEX  VARIABLE 

C  THIS  SUBROUTINE  IS  USED  TO  CALCULATE  THE  GI,HI  FUNCTION 

C 

COMPLEX  All, AI1P.BI1, BI1P 
COMPLEX  ZS , GI , GIP , HI ,HIP , ATI , AT2 
DIMENSION  AGT( 1001), AGP (1001) 

COMMON  ERR1 

PI-3.14159265369 

PI4-PI/4. 

SQPI-SQRT(PI) 

IARG-2 . 

IF (XS . GE . 4 . 7) IARG-3 
IF(IS.LE.-9.)IARG-1 
C 

C  IARG-1 :  NEGATIVE  LARGE  ARGUMENT;  THE  SOLUTION  IS  OSCILLATING 

C  IARG-2:  SMALL  ARGUMENT;  FOR  BOTH  POSITIVE  AND  NEGATIVE  ARGUMENT 

C  IARG-3:  POSITIVE  LARGE  ARGUMENT;  THE  SOLUTION  DAMPING  FAST 

C 
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Ti=ABS(XS) 

Tl-SQRT(Tl) 

XS14«SQRT(T1) 

XS32«T1*T1*T1 

ZS-CMPLI(XS,0.) 

CALL  AIBKZS, All, AliP.BIl, BI1P) 

GO  TO  ( 100 , 200 , 300 ) I ARC 
100  CONTINUE 

XS3-XS*XS*XS 

T2*-l . / (PI*XS)* ( 1 . +2 . /XS3+40 . / (XS3*XS3) ) 

T3-1 . / (PI*XS*XS) * (1 . +8 . /XS3+280 . / (XS3*IS3) ) 

HI-CMPLKT2.0.) 

HIP*CMPLX(T3,0.) 

GI-BI1-HI 

GIP-BI1P-HIP 

RETURN 

200  CONTINUE 

IF (ABS (IS) . NE . 0 . ) GO  TO  250 
GI-CMPLI(0 .  204976642478 ,  0 . ) 

GIP-CMPLX(0 . 149429462449,0. ) 

HI«2.*GI 
HIP-2 . *GIP 
RETURN 

260  CONTINUE 

C 

C  THE  INTEGRAL  IS  CALCULATED  USING  THE  DEFINITE  INTEGRAL  0  TO  AUG 

C  AND  THEN  ADD  UP  THE  REMAINDER  FROM  THE  REST  OF  THE  INTEGRAL 

C  USING  ASYMPTOTIC  EVALUATION 

C 

AUG-4. 

IF(XS.GE.3. )AUG-8. 

IF ( IS . LE . - 3 . ) AUG-2 . 

AUG2-AUG+AUG 

DELT-0.01 

NNT-AUG/DELT+1 
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DO  280  I«1,NNT 
XT*(I-1)*DELT 
ARGa-IT+IT*XT/3 ,+XS*XT 
T1*EXP(ARG) 

AGT(I)=EXP(ARG) 

AGP ( I ) *XT  *EXP ( ARG) 

280  CONTINUE 

Nl=(NNT-l)/2 

T2*0. 

T3=0. 

DO  290  J»1,N1 
J1=2*J 
J1M-J1-1 
J1P-J1+1 

T2=T2+(AGT(JlM)+4. *AGT(Jl)+AGT(JiP))*DELT/(3 . *PI) 
T3=T3+ (AGP (JIM) +4 . *AGP ( J 1 ) +AGP ( J1P ) ) *DELT/ (3 . *PI ) 

290  CONTINUE 

T4=XS*AUG-AUG*AUG2/3 . 

T6»EXP(T4) 

ISTT»1 . /(XS-AUG2) 

IST2-XSTT*XSTT 

XST4»IST2*XST2 

T7*T6*ISTT+ (-1 . +2 . +AUG*XST2- (2 . *XS+10 . *AUG2) *XST4) /PI 
T2«*T2+T7 

T8-T6*XST2* (1.-6. *AUG*XST2+ ( 8 . *IS+52 . *AUG2) *IST4) /PI 

T3*T3+4 . *T7+T8 

HI»CMPLI(T2,0.) 

HIP«CMPLX(T3,0.) 

GI-BI1-HI 

GIP-BI1P-HIP 

RETURN 

300  CONTINUE 

IS3»IS*IS+XS 

T2-1 . / (PI*IS) * (1 . +2 . /XS3+40 . / (IS3*XS3) ) 

T3—  1 .  / (PI*XS*IS) * ( 1 . +8 . /XS3+280 . / (XS3*IS3) ) 
GI»CMPLX(T2,0.) 

GIP"CMPLI(T3,0. ) 
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ARGa2.*XS32/3. 

T4aEXP(ARG) 

T5=T4/(SQPI*XS14) 

T6»T4*IS14/SQPI 

HI»CMPLI(T5,0.) 

HIP=CMPLX (T6 , 0 . ) 

RETURN 

END 


CCC— . - . - . . 

SUBROUTINE  FCTI(ID,FCT,X) 

COMPLEX  FII(8) ,FX(8) ,CJ ,FCT 
DIMENSION  XX (8) 

DATA  PI ,TPI ,SML/3 . 14159266 , 6.28318631,0.001/ 

DATA  IX/. 3.. 6,. 7.1. ,1.6, 2. 3, 4. ,6.6/ 

DATA  CJ/(0. ,1.}/ 

DATA  FX/ (0 . 6729 , 0 . 2677 ) , (0 . 6768 , 0 . 2682) , (0 . 7439 , 0 . 2649) , 

1(0. 8096 , 0 . 2322) , (0 . 873 , 0 . 1982) , (0 . 9240 , 0 . 1677 ) , (0 . 9668 ,0.1073), 

2(0.9797,0.0828)/ 

DATA  FIX/ (0 . , 0 . ) , (0 . 6195 , 0 . 0025) , (0 . 3366 , -0 . 0666) , 
1(0.2187,-0. 0767) .(0.127,-0. 068) , (0 . 0638 , -0 . 0606) , 

2(0. 0246, -0.0296), (0.0093, -0.0163)/ 

IF(X.GT.6.6)G0  TO  1 
IF (I. GT. 0.3) GO  TO  10 
C 

C  ID-1  DIFFRACTION  COEFFICIENT  FX 

C  ID-2  SLOPE  DIFFRACTION  COEFFICIENT  FIS 

C 

C!!!  SMALL  ARGUMENT  FORM 

FCT- ( (1 . 253 , 1 . 263) *SQRT (X) - (0 . ,2.)*X-0.6667*X*I)*CEXP(CJ*I) 

IF(ID.EQ.2)  FCT-2 . *CJ*X* ( 1 • -FCT) 

RETURN 

CH!  LINEAR  INTERPOLATION  REGION 

10  DO  11  N-2,7 

11  IF(I.LT.IX(N))GO  TO  12 
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12  FCT=FXX(N)*(X-XI(N))+FX(N) 

IFCID.EQ.2)  FCT=2 . ( 1 . -FCT) 

RETURN 

C! ! !  LARGE  ARGUMENT  FORM 

1  IF(ID.EQ.l)  FCT=1. +CMPLX( -0.75/1,0. 5) /X 

IF < ID . EQ .  2)  FCT*1 .  +CMPLX  (-3 .  76/X ,  1 . 5)  /?. 
RETURN 
END 
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FUNCTION  COMPLETE  AIRY 

These  functions  compute  values  for  the  Complete  Airy  functions. 

complex  function  ai  (  z  ) 

complex  z,  aip,  bi,  bip,  zl,  z2,  z3,  z4 

call  aibi  (  z,  zl,  z2,  z3,  z4  ) 

ai  =  zl 

return 

entry  aip 

call  aibi  (  z,  zl,  z2,  z3,  z4  ) 

aip  *  z2 

return 

entry  bi 

call  aibi  (  z,  zl,  z2,  z3,  z4  ) 

bi  *  z3 

return 

entry  bip 

call  aibi  (  z,  zl,  z2,  z3,  z4  ) 

bip  *  z4 

return 

end 

SUBROUTINE  AIBI(Z,AI,AIP,BI,BIP) 

COMPLEX  Z,AI,AIP,B1,BIP 
IF (CABS (Z) .GT.6 . )G0  TO  12 
CALL  AIBI1(Z,AI,AIP,BI,BIP) 

RETURN 

12  CALL  AIBI2(Z,AI,AIP,BI,BIP) 

RETURN 

END 

SUBROUTINE  AIBI1(Z,AI,AIP,BI,BIP) 

C  THIS  PROGRAM  CALCULATES  THE  AIRY  FUNCTIONS  AI(Z),BI(Z), 
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C  AND  THEIR  DERIVATIVES  AIP(Z) ,BIP(Z) , 

C  REF.  ABRAMOWITZ  AND  STEGUN,  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS. 

C  FOR  CABS(Z)  .LE.  6.0  ,A  TAYLOR  SERIES  IS  USED. 

C  ARG(Z)  MAY  TAKE  ANY  VALUE.  SEE  (10.4.2)  TO  (10.4.6)  . 

COMPLEX  Z,AI,AIP,BI,BIP 
COMPLEX+16  F.G.FP.GP 
DOUBLE  PRECISION  CC1.CC2 

DATA  S3, CC1.CC2/1. 732050808,. 365028063887817, .268819403792807  / 

CALL  FZGZ(Z,F,G,FP,GP) 

AI=CC1*F-CC2*G 

AIP«CC1*FP-CC2*GP 

BI»S3+(CC1*F+CC2*G) 

BIP-S3* (CC1*FP+CC2*GP) 

RETURN 

END 

SUBROUTINE  FZGZ(Z,F,G,FP,GP) 

C  THE  AUXILIARY  FUNCTIONS  F(Z) ,G(Z) ,FP(Z) ,GP(Z)  ARE  COMPUTED  AS  IN 

C  "TABLES  OF  THE  MODIFIED  HANKEL  FUNCTIONS  OF  ORDER  ONE-THIRD  AND 

C  OF  THEIR  DERIVATIVES", COMPUTATION  LAB,  HARVARD  UNIV.  PRESS, 1946. 

C0MPLEX*16  F,G,FP ,GP,Z3,Z3M,ZD 
COMPLEX  Z 

REAL *8  AM , BM , CM , DM , AO , BO , CO , DO 
REAL  ZMBD(6) 

INTEGER  MAX (6) 

DATA  ZMBD  /6.1,  6.6,  4.8,  4.1,  3.2  / 

DATA  MAX  122,  19,  16,  14,  11  / 

ZD-O. DO 
ZD-Z 
AO-1. DO 
BO-1. DO 
C0-O.6DO 
DO-1. DO 
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Z3=(ZD**3)/200 

Z3K=Z3 

ZMAG*CABS(Z) 

DO  3  M*1 ,6 

3  IFCZMAG  .LE.  ZMBD(M) )MADMAX=MAX(M) 

FxAO 

G=BO 

FP*CO 

GP=DO 

DO  10  M=1 .MADMAI 
TM«FL0AT(3*M) 

AMS1200 .  DO+AO/TM/  (TM-1) 

BM*200 . DO*BO/TM/ (TM+1) 

CM»200 . DO+CO/TM/ (TM+2) 

DM=200 .DO+DO/TM/ (TM-2) 

F=F+AM*Z3M 

G'=G+BM*Z3H 

FP»FP+CM+Z3M 

GP*=GP+DM*Z3M 

Z3M=Z3M*Z3 

AO=AM 

BO*BH 

CO®CM 

DO*DM 

10  CONTINUE 

G»ZD*G 

FPx(ZD**2)*FP 

RETURN 

END 

SUBROUTINE  AIBI2(II,AI,AIP,BI,BIP) 

C  THIS  PROGRAM  CALCULATES  THE  AIRY  FUNCTIONS  AI(II) ,BI(XI) , 

C  AND  THEIR  DERIVATIVES  AIP(XX) ,BIP(II) . 

C  REF.  ABRAMOWITZ  AND  STEGUN,  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS. 

COMPLEX  Z,AI,AIP,BI,BIP,XX 
COMPLEX  Z2S , ZTB , ZT , ZT2 , ZT3 , ZT4 , ZT5 
COMPLEX  CT1,A2L2,EIPI3.EIPI6,C,S 
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DATA  RTPI , TWORPI , RTOP , POF 

%  / 1.772453861, 3. 54490T702, .797884661 , . 785398164  / 

DATA  A2L2 , EIPI6 , EIP 13 

y,  /(O... 346573590), (.866025404,. 5), (.5,. 866026404)  / 
DATA  Cl/ . 069444444/ , C2 / . 037133487/ , C3/ . 037993059/ , 

%  C4/. 057649190/, C5/. 116099064/ 

DATA  Dl/- . 097222222/ , D2/- . 043885030/ , D3/- . 042462830/ , 
*/.  D4/- .  062662163/ ,  D5/- .  124105896/ 

ZTB«(2 . /3 . ) +11**1 . 6 

ARG3 ATAN2 ( AIMA6 (XI ) .REAL (IX)) 

IF(ABS(ARG) .GE.2.1)  GO  TO  100 

C  EQN.  (10. 4. 59), (10. 4. 61) 

Z26=IX+*.25 

ZT*ZTB 

ZT2**ZT*ZT 

ZT3=ZT2*ZT 

ZT4=ZT2*ZT2 

ZT6«ZT3*ZT2 

CT1-CEIP ( -ZT) /TWORPI 

AI«CT1/Z25*(1-C1/ZT+C2/ZT2-C3/ZT3+C4/ZT4-CS/ZT6) 

AIP— CT1*Z25*(1-D1/ZT+D2/ZT2-D3/ZT3+D4/ZT4-D6/ZT5) 
IF(ARG .LT.O. )G0  TO  20 
ZT«(0. ,-l . )*ZTB 

C  EQN.  (10. 4. 65), (10. 4. 68)  WITH  UPPER  SIGNS. 

Z*XX/EIPI3 
CT1-ZT+P0F-A2L2 
BI-EIPI6 
BIP-1./EIPI6 
GO  TO  30 

20  ZT*(0.,1.)*ZTB 

C  EQN.  (10. 4. 65), (10. 4. 68)  WITH  LOWER  SIGNS. 

Z»IX*EIPI3 

CT1«ZT+P0F+A2L2 

EI-1./EIPI6 

BIP*EIPI8 

30  S-CSIN(CTl) 

C«CC0S(CT1) 
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Z25=Z** .26 

ZT2*ZT*ZT 

ZT3=ZT2*ZT 

ZT4=ZT2*ZT2 

ZT6*ZT3*ZT2 

BI*BI*RT0P/Z25*(S+(l-C2/ZT2+C4/ZT4)-O(Cl/ZT-C3/ZT3+C6/ZT6)) 
BIP*BIP*RT0P+Z25*(O(l-D2/ZT2+D4/ZT4)+S*(Dl/ZT-D3/ZT3+D6/ZT6) ) 
RETURN 

100  ZT=(0. , l.)*ZT8 

C  EQN.  (10.4.60) , (10.4.62) , (10.4.64) ,(10.4.67) 

IF ( ARG . LT . 0 . ) ZT*-ZT 
Z=-XX 

Z26*Z**.26 

ZT2°ZT+ZT 

ZT3«ZT2*ZT 

ZT4-ZT2+ZT2 

ZT6°ZT3*ZT2 

CT1-ZT+P0F 

S-CSIN(CTl) 

OCCOS(CTi) 

11*1 . /RTPI/Z25* (S* (1-C2/ZT2+C4/ZT4) -C* (C1/ZT-C3/ZT3+C6/ZT5) ) 

AIP*-Z26/RTPI*(0(1-D2/ZT2+D4/ZT4)+S+(D1/ZT-D3/ZT3+D6/ZT6)) 

BI»1 ./RTPI/Z26*(C*(1-C2/ZT2+C4/ZT4)+S*(C1/ZT-C3/ZT3+C6/ZT6) ) 

BIP-Z26/RTPI* (S* ( 1-D2/ZT2+D4/ZT4) -C* (D1/ZT-D3/ZT3+D6/ZT5) ) 

RETURN 

END 
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